arXiv:1501.06094v2 [math.PR] 15 Sep 2015 


Wavelet estimation for operator fractional Brownian motion 


Patrice Abry Gustavo Didier 

Physics Lab Mathematics Department 

CNRS and Ecole Normale Superieure de Lyon Tulane University 

September 17, 2015 


Abstract 

Operator fractional Brownian motion (OFBM) is the natural vector-valued extension of 
the univariate fractional Brownian motion. Instead of a scalar parameter, the law of an OFBM 
scales according to a Hurst matrix that affects every component of the process. In this paper, 
we develop the wavelet analysis of OFBM, as well as a new estimator for the Hurst matrix 
of bivariate OFBM. For OFBM, the univariate-inspired approach of analyzing the entry-wise 
behavior of the wavelet spectrum as a function of the (wavelet) scales is fraught with difficulties 
stemming from mixtures of power laws. The proposed approach consists of considering the 
evolution along scales of the eigenstructure of the wavelet spectrum. This is shown to yield 
consistent and asymptotically normal estimators of the Hurst eigenvalues, and also of the 
coordinate system itself under assumptions. A simulation study is included to demonstrate 
the good performance of the estimators under finite sample sizes. 


1 Introduction 

An M"^-valued stochastic process {A(t)}tg]R is said to be operator self-similar (o.s.s.) when its law 
scales according to a matrix (Hurst) exponent H, i.e., 

{X{ct)}teR = {c^X{t)}teR, c> 0, (1.1) 

where {c)H^/kl and = denotes the equality of finite-dimensional distributions. 

No specific assumption on the eigenstructure of H is imposed, e.g., canonical vectors are not 
necessarily eigenvectors. The notion of operator self-similarity underpins the natural multivariate 
generalization of the univariate fractional Brownian motion (FBM): an operator fractional 
Brownian motion (OFBM) Bh = {BHit)}t£R is a proper Gaussian, o.s.s., stationary increment 
stochastic process. In this paper, we propose using the wavelet eigenstructure of OFBM to 
estimate H. The main motivation behind this methodology is to avoid the difficulties stemming 
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from the extrapolation of univariate techniques to a multivariate context, where operator scaling 
laws such as (1.1) may arise. 

Inferential theory for univariate self-similar processes now comprises a voluminous and well- 
established literature. A non-exhaustive list includes Fox and Taqqu (1986) and Robinson (1995o, 
19956) on Fourier domain methods, and Wornell and Oppenheim (1992), Flandrin (1992), and 
Veitch and Abry (1999) on wavelet domain methods, among many others. The multivariate 
framework evokes several applications where matrix-based scaling laws are expected to appear, 
such as in long range dependent time series (Marinucci and Robinson (2000), Davidson and 
de Jong (2000), Chung (2002), Dolado and Marmol (2004), Davidson and Hashimzade (2008), 
Kechagias and Pipiras (2015)) and queueing systems (Konstantopoulos and Lin (1996), Majewski 
(2003, 2005), Delgado (2007)). Like FBM in the univariate setting, OFBM is a natural starting 
point in the construction of estimators for operator self-similar processes due to its tight connection 
to stationary fractional processes and its being Gaussian (on the general theory of o.s.s. processes, 
see Laha and Rohatgi (1981), Hudson and Mason (1982), Maejima and Mason (1994), Cohen et 
al (2010)). 

A characterization of the covariance structure of OFBM can be derived from stochastic integral 
representations. Under a mild condition on the eigenvalues of the exponent H (see (2.9)), Didier 
and Pipiras (2011) showed that any OFBM Bh admits a harmonizable representation 

= I [ + xZ^Widx)} (1.2) 

1 Jr ZX ) tSK 

for some complex-valued matrix A. In (1.2), x± = max{±x,0}, 

D = H-^I, (1.3) 

and B{dx) is a complex-valued random measure such that B{—dx) = B{dx), EB{dx)B{dx)* = dx, 
where * represents Hermitian transposition. Expression (1.2) shows that the law of an OFBM 
can be fully described based on the scaling matrix H and the spectral parameter A (see Remark 
2.1 on the parametrization). 

Let H = PJ}jP~^, P G GL(n,C), be the Jordan form of the Hurst parameter in (1.2) 
(see Section 2 for matrix notation). If H is diagonal, then we can assume that P takes the 
form of a scalar matrix P = pi, where p G M and I is the n x n identity matrix, and that 
Jh = diag(/ii,..., hn)- In this case, (1.1) breaks down into simultaneous entry-wise expressions 

{X{ct)}tm = {{c^^X^{t),...,c’^-Xn{t)y}teR, C>0. (1.4) 

Relation (1.4) is henceforth called entry-wise scaling. In particular, under (1.4) an OFBM is 
a vector of correlated FBM entries (Amblard et al. (2012), Coeurjolly et al. (2013)). Several 
estimators have been developed by building upon the univariate, entry-wise scaling laws, e.g., 
the Fourier-based multivariate local Whittle (e.g., Shimotsu (2007), Nielsen (2011)) and the mul¬ 
tivariate wavelet regression (Wendt et al. (2009), Amblard and Coeurjolly (2011), Achard and 
Gannaz (2014)). However, if H is non-diagonal, i.e., if P is not a scalar matrix, then the relation 
(1.1) mixes together the several entries of X. The estimation problem under a non-scalar P turns 
out to be rather intricate and calls for the construction of methods that are multivariate from 
their inception. 

Although the emergence of o.s.s. processes in applications is rightly expected - e.g., as func¬ 
tional weak limits of multivariate time series -, there is no specific reason to believe a priori that 
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scaling laws occur predominantly entry-wise and exactly along the canonical axes. Indeed, this 
is palpably not true in several applications such as fractional blind-source separation (see Didier 
et al. (2015)) and fractional cointegration (see Robinson (2008) for a bivariate local Whittle es¬ 
timator). The framework of multivariate mixed fractional time series subsumes both cases. Let 
W = {Wt}t£Z be an unobserved signal. To fix ideas suppose that the signal is an operator frac¬ 
tional Gaussian noise (OFGN), namely, it has the form Wt = — 1) with spectral 

parameter A\y (see (1.2)). Further assume that H\y = diag(/ii,..., hn)- The observed signal has 
the form Y = {Yt}t£Z = {PWt}t&z, for a non-scalar, mixing matrix P £ GL(n,M). Then, the 
mixed process Y is another OFGN with parameters Hy = Pdiag(/ii,..., hn)P~^ and Ay = PAw- 
In blind source separation, the entries of W are uncorrelated, whereas in cointegration they are 
typically correlated. 

From a mathematical standpoint, the mixing of scaling laws can be illustrated by means of the 
expression for the spectral density fx of an OFBM with Hurst parameter H = i-’diag(/ii, h 2 )P~^-, 
0 < hi < /i 2 < 1, R G GL(2,M) (see condition (2.12)). For M := P~^AA*{P*)~^ £ M{n) (c.f. 
condition (2.11)) and x > 0, the spectral density takes the form {jx{x)ij^ = X ^AA*x , 
where 

fx{x)n =phmnx~^‘^^ + 2pnPi2mi2X~^'^^^'^^^ + 

fx{x)21 = PllP2imuX~‘^'^^ + {puP22 + Pl2P2l)mi2X~^‘^^+^‘^'> + Pl2P22m22X~‘^'^^ , 

fx{x)22 = -f 2p2lP22mi2X“("'l+*) -f P22"i22X"^'^^ 

and di = hi — 1/2, d 2 = h 2 — 1/2 (see (4.18) for the analogous expression in the wavelet domain). 
The univariate-inspired approach of setting up a Fourier-domain log-regression - e.g.. Whittle-type 
estimators - has to cope with the double-sided challenge of mixed power laws. On one hand, under 
mild assumptions on the amplitude coefficients, the dominant power law always prevails 

around the origin of the spectrum. On the other hand, and paradoxically, even if the estimation 
of d2 is the target, the magnitude of the amplitude coefficients themselves can arbitrarily bias the 
estimate over finite samples by masking the dominant power law. 

In this work, we build upon the wavelet analysis of (n-dimensional) OFBM to propose a 
novel wavelet-based estimation method for bivariate, and potentially multivariate, OFBM. The 
method yields the Hurst eigenvalues of H and, under mild assumptions, also its eigenvectors when 
P £ 0(2). Its essential ingredient, and the main theme of this paper, is a change of perspective; 
instead of considering the entry-wise behavior of the wavelet spectrum as a function of wavelet 
scales, it draws upon the evolution along scales of the eigenstructure of the wavelet spectrum. 
This way, it avoids much of the difficulty associated with inference in the presence of mixed power 
laws, as we now explain. 

For a wavelet function ip £ L^(M) with a number N-^ of vanishing moments (see (2.6)), the 
(normalized) vector wavelet transforms of OFBM is naturally defined as 

9 D{2^, k) = 2-P^ [ 2-P^'iP{2-H - k)BH{t)dt, j G N U {0}, A: G Z, (1.5) 
Jr 

provided the integral in (1.5) exists in an appropriate sense. The wavelet-domain process 
{D{2^, k)}k£Z is stationary in k and o.s.s. in j (Proposition 3.1). Moreover, whereas the orig¬ 
inal stochastic process Bnit) displayed fractional memory, the covariance between (multivariate) 
wavelet coefficients decays as a function of \2^k — 2^ k'\ according to an inverse fractional power 
controlled by N-^ (Proposition 3.2). The wavelet spectrum (variance) at scale j is the positive 
definite matrix 

ED{2^,k)D{2^,ky = ED{2yO)D{2yoy =: RIF(2^), 
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and its natural estimator, the sample wavelet variance, is the matrix statistic 


W{2^) 


K, 


Y,Di2^,k)D{2\kr, 

k=l 


~ 2i ’ j ~ ■ ■ ■ 1 jrn, 


( 1 . 6 ) 


for a dyadic total of i/ (wavelet) data points. Within the bivariate framework, the univariate-like 
entry-wise scaling approach would consist of exploiting the behavior of each component W ( 2 '^)j^^j 2 ) 
ii,Z 2 = 1,2, of the sample wavelet transform W{2^) as a function of the scales 2^. Apart from 
an amplitude effect, the entries are then controlled by the dominant Hurst eigenvalue h 2 (see 
expression (4.16)). Figure 1, top panels, illustrates the fact that this precludes the estimation of 
hi. 

The proposed estimators of the Hurst eigenvalues hi and h 2 are 


^i(2^) 


logAi(2-^) 
21og(2i) ’ 


^ 2 ( 2 ^') 


logA2(2^) 

21og(2i) ’ 


(1.7) 


where Ai(2-^) < A2(2-^) are the eigenvalues of the positive definite symmetric matrix W{2^) (see 
Definition 4.1 for the precise assumptions). However, as usual with operator self-similarity, the 
finite sample expressions for Ai(2'^) and A 2 ( 2 -^) themselves involve a mixture of distinct power 
laws <; /J 2 ). For this reason, one must take the limit at coarse scales, 

namely, the scale itself must go to infinity. It is a remarkable fact that the power law ends 
up prevailing in the expression for Ai(2'^) (see Figure 1, bottom panels, and the striking contrast 
with the top panels; see Remark 4.1 for a mathematically motivated, intuitive discussion). The 
convergence of (1.7) in turn allows for the convergence of associated sequences of eigenvectors when 
P is orthogonal. Moreover, simulation studies show that the estimation procedure is accurate and 
computationally fast. The asymptotics are mathematically developed in two stages. In the first, 
the wavelet scales (octaves) are held fixed and the asymptotic distribution of the sample wavelet 
transform is obtained (Proposition 3.3 and Theorem 3.1). In the second, one takes the limit with 
respect to the scales themselves. However, the latter must go to infinity slower than the sample 
size, a feature that our estimators share with Fourier or wavelet-based semiparametric estimators 
in general (e.g., Robinson (1995a), Moulines et al. (20076, 2007a, 2008)). 

Our results are related to the literature on the estimation of operator stable laws via eigenvalues 
and eigenvectors of sample quadratic forms (see Meerschaert and Scheffler (1999, 2003)). In this 
context, one encounters the same problem with the prevalence of some dominant power law (i.e., 
the tail exponent) in most directions. In Becker-Kern and Pap (2008), a similar philosophy is 
applied in the time domain to produce one of the very few available estimators for authentic, mixed 
scaling o.s.s. processes of dimension up to 4. However, the asymptotics provided are restricted to 
consistency. In our work, the wavelet transform is the main tool for ensuring the consistency and 
asymptotic normality of the proposed estimators. 

The paper is organized as follows. Section 2 contains the notation, assumptions and basic 
concepts. Section 3 is dedicated to the wavelet analysis of n-dimensional OFBM, as well as the 
asymptotics of the wavelet transform for fixed scales (most of the proofs can be found in Section 
B). In Section 4, the estimation method for the Hurst exponent of bivariate OFBM is laid out 
in full detail and its asymptotics are established at coarse scales. Section 5 displays finite sample 
computational studies, including one of the performance of the estimators under blind source 
separation and cointegrated instances, with the purpose of illustrating the robustness of the 
estimators’ performance with respect to different parametric scenarios. The research contained 
in this paper leads to a number of interesting open questions, which are mentioned in Section 6. 
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Among these is the extension of the consistency and asymptotic normality to any dimension n > 3, 
which will generally require dispensing with explicit formulas for eigenvalues and eigenvectors. 
The appendix contains several auxiliary mathematical results. In addition, in Section C, the 
performance of the estimators is established under the assumption that only discrete observations 
are available, instead of a full sample path as in (1.5). 



j = log2 2' 





Figure 1: Entry-wise vs eigenvalue-based estimation. From one synthetic realization of 
OFBM, the top row displays (black solid lines with =t=), from left to right, log 2 kF(2'^)i^i vs j, 
log 2 W vs j and log 2 W{ 2 ^) 2,2 vs. j, on which the asymptotic behavior j 2 h 2 is superimposed 
(red dashed line with ‘ 0 ’). All auto- and cross-components are then driven by the dominant Hurst 
eigenvalue h 2 , which precludes the estimation of the eigenvalue hi. The first two bottom row 
plots display (black solid lines with *), from left to right, log 2 Ai( 2 -^) vs j and log 2 A 2 ( 2 -^) vs j, 
with their respective asymptotic trends j2hi and j 2 h 2 superimposed (red dashed line with ‘ 0 ’). 
This demonstrates that both Hurst eigenvalues hi and h 2 can be estimated (see Section 5.1 for 
the simulation details). The bottom right plot displays the increments of the generated OFBM. 


2 Notation and assumptions 

All through the paper, the dimension of OFBM is denoted by n > 2. 

We shall use throughout the paper the following notation for finite-dimensional operators 
(matrices). All with respect to the field M, M(n) or M(n,M) is the vector space of all n x n 
matrices (endomorphisms), GL{n) or GL(n,M) is the general linear group (invertible matrices, 
or automorphisms), 0{n) is the (orthogonal) group of matrices O such that 00* = I = 0*0 
(i.e., the adjoint operator is the inverse), and S{n, M) is the space of symmetric matrices. We also 
write In to indicate the dimension of the identity matrix I. A block-diagonal matrix with main 
diagonal blocks Vi,... ,Vn oi m times repeated diagonal block V is represented by 

diag('Pi,...,Pn), diag^(P), (2.1) 
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respectively. The symbol || • || represents a generic matrix or vector norm. For q G N, the 
entry-wise norm of an m x n real-valued matrix M is denoted by 

m n , 

= ' (2.2) 
il=l 12=1 

and ||M||;oo = supj^ ^ generic matrix M G M{n, C) has real and imaginary parts 3?(M) 

and 3'(M), respectively. The functions 

'^i{v),7ri^,i2{M), V G M"-, M G M{n,R), (2.3) 

denote, respectively, the i-th projection (entry) of the vector v and the (zi,i 2 )-th projection 
(entry) of the matrix M, ii,i 2 = l,...,n. For S = £ 5(n,M), we define the 

operator 

^ 605 ( 5 ") — (sil, S 12 , . . . , Sin; 'S22) • • • ) S2n: • • • > -Sn—l,n—1; Sn—l,nj ■ (^•^) 

In other words, vecs{-) vectorizes the upper triangular entries of the symmetric matrix S. 

When establishing bounds, C stands for a positive constant whose value can change from one 
line to the next. For a sequence of random vectors {X;, P{Yi = 0) = 0, we write 

Xi ~ Yi (2.5) 

p 

to mean that XijYi —)■ 1, ? —)■ oo. Note that this does not imply that converges in 

probability. Relations of the type (2.5) will often appear in the proofs of the results in Section 4. 
We write X = Y when the random vectors X and Y have the same distribution. 

All through the paper, we will make the following assumptions on the underlying wavelet 
basis. For this reason, such assumptions will be omitted in the statements. 

Assumption (VFl): ip G L^(M) is a wavelet function, namely. 


ip‘^{t)dt = l, / Pip{t)dt = 0, q = 0,1,..., — 1, N^>2. 


Assumption (VF2): 
Assumption (VF3): there is 


supp('0) is a compact interval. 
a > 1 such that 

sup \ip{x)\{l -I- |x|)“ < oo. 


( 2 . 6 ) 


(2.7) 


( 2 . 8 ) 


Under (2.6), (2.7) and (2.8), ip is continuous, ip{x) is everywhere differentiable and its first 
— 1 derivatives are zero at x = 0 (see Mallat (1999), Theorem 6.1 and the proof of Theorem 
7.4). 

Example 2.1 If i/; is a Daubechies wavelet with vanishing moments, supp('!/;) = [0, 2N^ — 1] 
(see Mallat (1999), Proposition 7.4). 

Starting from the harmonizable representation (1.2), throughout the paper we will make the 
following assumptions on the OFBM Bh = 
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Assumption (OFBMl): the eigenvalues (characteristic roots) of the matrix exponent H 
satisfy 


0<»(/ifc)<l, A: = l,...,n. (2.9) 

Assumption (OFBM2): 

iR(AA*) has full rank. (2-10) 

The condition (2.9) generalizes the familiar constraint 0 < < 1 on the Hurst parameter of 

a FBM. As shown in Didier and Pipiras (2011), it ensures the existence of the harmonizable 
representation (1.2). Also, (2.9) implies that the OFBM under consideration has mean zero, 
which follows by the same reasoning as in Taqqu (2003), p.7, property {ii). In turn, recall that 
a stochastic process is called proper when the distribution of X{t) is full dimensional for t ^ 0. 
The condition (2.10) is sufficient (though not necessary) for the integral on the right-hand side of 
(1.2) to be a proper stochastic process and hence to define an OFBM. 

The next two assumptions will appear in some of the results. 

Assumption (OFBM3): 

3=(AA*) = 0. (2.11) 

Assumption (OFBM4): Bh = is a bivariate OFBM with scaling matrix 

H = PJhP~^ = Pd\a,g{hi,h2)P~^■: 0 </ii </12 < 1, PgGL(2,M), (2-12) 

where the columns of P are unit vectors. 

The condition (2.11) is equivalent to time reversibility, namely, {Bni—t)} = {Bnit)}. In turn, 
the latter is equivalent to the existence of a closed form expression for the covariance function, 
i.e., 

EBH{s)BH{t)* = ^{\t\^E\t\^* + -\t- s\^E\t - s|^*}, (2.13) 

where S := FiH/f(l)i?//(l)* (Didier and Pipiras (2011), Proposition 5.2). Time reversibility is 
used in some of the results in Section 3 and in all Section 4. The bivariate framework (2.12) is 
only used in limits at coarse scales (Section 4), due to the availability of convenient formulas for 
eigenvalues and eigenvectors. 

Remark 2.1 In regard to the parametrization, the connection between S and {P[,AA*) can be 
obtained implicitly based on the second moment of the harmonizable representation (1.2) at t = 1, 
and it can be worked out explicitly under (2.13) and stronger additional conditions, e.g., assuming 
H is diagonalizable with real eigenvalues and eigenvectors. 

The condition (2.12) renders OFBM identifiable, namely, the mapping from the parametriza¬ 
tion H into the space of OFBM laws is injective for a fixed spectral parameter AA*. This is a 
consequence of only considering Hurst matrices H with real eigenvalues. For a general discussion 
of the (non)identifiability of OFBM, see Didier and Pipiras (2012). 

3 Wavelet analysis 

In this section, we carry out the wavelet analysis of n-dimensional OFBM. The proof of Propo¬ 
sition 3.1 can be found in Section A, whereas Section C contains the proofs of Proposition 3.2, 
Proposition 3.3 and Theorem 3.1. 
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3.1 Basic properties 

The normalized wavelet transform (1.5) is itself a vector-valued random field in the scale and 
shift parameters j and k, respectively. It will be convenient to make the change of variables 
z = 2~H — k, and reexpress 

D{2\k)= [ 'ip{z)BH{2^z + 2^k)dz. (3.1) 

Jr 

As in the univariate case, the wavelet coefficients of OFBM exhibit a number of nice properties. 
The next proposition describes such properties as well as the general form of the wavelet spectrum 
(variance). 

Proposition 3.1 Under the assumptions (0FBMl)-(2), let {D{2^,k)}o,s in (3.1). 
Then, 

(PI) the wavelet transform (1.5) is well-defined in the mean square sense, and ED{2fik) = 0; 
(P2) (stationarity for a fixed scale) {D{2fi k + /i)}fcez = {D{2fi k)}kez, h G Z; 

(P3) (operator self-similarity over different scales) {D{2fi k)}kez = D{l,k)}kezi 
(P4) for some C > 0, the wavelet spectrum EW{2fik) = EW{2h) is given by 

EW{2^) = C [ (3.2) 

Jr X 


(P5) the wavelet spectrum satisfies the operator scaling relation 

EW{2^) = 2^^EW{1)2^^* = 2^^{ f f fi{t)fi{t') EBH{t)BH{t'ydtdt'\2^^\ 

^ Jr Jr ’ 

j G N. In partieular, under (2.11), 

EW{2y = -I- 2^^{ f f fi{t)fi{t') \t - t'\^E\t - t'\^*dtdt'\2^^* ; (3.3) 

2 Wr Jr 1 


(P6) in analogy to (P5), 1T(2-^) = 2^^{^ k)*)2^^*, where Kj is given in (1.6); 

(P7) the wavelet spectrum has full rank, namely, det EW{2^) 0, j £ N. 

Fix some 0 < <5 < 1/2, and consider the range of wavelet parameters j, f, k, k' such that 

max{2-^ , 2-^^} 6 

—^—- < -. 13.4i 

2l/c — 2Tk' max{l, length(supp(V’))} 

If the parameters (j, k) and {j', k') of two wavelet coefficients satisfy (3.4), then we can interpret 
that the latter are “far apart” in the parameter space. The next result provides a notion of decay 
of the covariance between wavelet coefficients under (3.4). The proof is similar to that for the 
univariate case, but we provide it in Section B for the reader’s convenience. 



Proposition 3.2 Under the assumptions (OFBMl)-(3) and (3.4), the eovarianee between 
wavelet coeffieients (3.1) satisfies the relation 

ED{2fik)D{2^',k')* 


-2^'k'\^(0^X(^)) \2^k-2^'k'\^*}, (3.5) 


where ( 


'2\2ik-2i'k'\‘^^^ 

is an entry-wise bounded symmetrie-matrix-valued function that de- 

ii ,i2=l,...,n 


pends only on N^. As a consequence, 

\\ED{2fik)D{2^',kr\\<C{N^,j,j') 


\log'^ \ 2^ k - 2^' E\ 


\2^k — 5R(^) ’ 


(3.6) 


where k is the dimension of the largest Jordan block in the spectrum of H, and C{N,p, j, j') = 
C{N^) for some constant C{N^) > 0. 


3.2 Asymptotics for sample wavelet transforms: fixed scales 

As typical in the asymptotic study of averages, we begin by investigating the asymptotic covariance 
of the sample wavelet transforms W{2^). 

For FBM, the asymptotic covariance between wavelet transforms 1F(2-^), VF(2'^ ) G M is not 
available in closed form since it depends on the wavelet function, which is itself not available in 
closed form (c.f. Bardet (2002), Proposition II.3). Operator self-similarity adds a layer of intricacy, 
since in general exact entry-wise scaling relations are not present. 

For notational simplicity, let 

X* = (Ai,...,A0 = idi{2fik),...,dn{2fik)), Y* = iYu...,Yn) = idi{2^\k'),... ,dn{2^',k')), 

(3.7) 

where di{2fik), i = l,...,n, is the i-th entry of the wavelet transform vector D[2fih). The 
bivariate case serves to illustrate the computation of covariances. 

Example 3.1 For a zero mean, Gaussian random vector Z G M™, the Isserlis theorem (e.g., 
Vignat (2012)) yields 

E{Zi... Z2k) = Efi E{ZiZj), E(Zi ... Z 2 fc+i) = 0, k = 1,... ,[m/2\. (3.8) 


The notation ^ ](([ stands for adding over all possible /c-fold products of pairs E{ZiZj), where the 
indices partition the set 1,..., 2k. So, let X and Y be as in (3.7) with n = 2. Then, 



/ 2{E{XiYi)f 2{E{XiY2)f 2E{XiY{)E{XiY2) \ 

2{E{X2Yi)f 2{E[X2Y2)f 2E{X2Yi)E{X2Y2) , (3.9) 

V 2E(AiYi)E(A2Yi) 2E(Aiy2)^(X2y2) C33 / 

where C33 = E{XiYi)E{X2Y2) + E{XiY2)E{X2Yi). 
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Expression (3.9) shows that the asymptotic behavior of the second moments of W{2^) involves 
several cross products. A notationally economical way of tackling this difficulty is by resorting to 
Kronecker products. For instance, in the bivariate case, 


M(4, M) 9 EXY* ® EXY* 


( E{XiYi) E{XiY2) \ 

V E{X2Yi) E{X2Y2) ) 


( E{XiYi) E{XiY2) \ 

V E{X2Yi) E{X2Y2) ) 


contains all the 9 terms (two-fold products of cross moments), as well as a few repeated ones, 
needed to express (3.9). In view of (3.8), this fact extends to general dimension n by means of 
the relations 


CoY{Xi,Xi,,Y,,Y,,) = E{X,,Yj,)E{Xi,Y,,) + E{Xi,Y,,)E{Xi,Yj,), (3.10) 

for = l,...,n. The next proposition provides an expression that encompasses the 

asymptotic fourth moments of the wavelet coefficients. 

Proposition 3.3 Let Bh = {BH{t)}t£R be an OFBM under the assumptions (OFBMl) - (3) 
and let Kj be as in (1.6). As v ^ oo, for every pair of octaves j, j', 

{i) 

Kj K-, 

ED{2\ k)D{2^', k')* ® ED{2\k)D{2^\k')* 

^ ^ k=l k'=l 

OO 

2/) j: Lljj/{z) 0 Lljji{z), (3-11) 

z=—oo 

where 

Q,jj/{z) :=— f f ^p{t)^p{t')\{2H—2^ t')+z gcd{2^,2^ )\^'E\{2H—2^ t')+z gcd{2^,2^ )\^ dt'dt 

(ii) there is a matrix Gjj> G M{n{n + 1)/2,M), not necessarily symmetric, such that 

Cov(vec5lT(2^'),vec5lT(2^')) ^ Gjf, 

where the entries ofGjj> can be retrieved from (3.11) by means of (3.10) (see (2.4) on the 
notation vecs). 

Remark 3.1 The definition of wavelet only requires > 1, but > 2 (see (2.6)) is needed 
for the convergence in Proposition 3.3 (see (B.32)). 

The next theorem establishes the asymptotics for the vectorized sample wavelet transforms 
(vec 5 lT( 2 -^))j=j^^..,j^ at a fixed set of octaves. 

Theorem 3.1 Let Bh = {BH}t£R be an OFBM under the assumptions (OFBMl)-(3). Let 
ji < ... < jm be a fixed set of octaves. Then, 

(v^(vec5lT(2^)-vec5AIT(2^'))) , 4A^^^(0,F), (3.12) 

as n ^ oo (see (2.4) on the notation vecsj- In (3.12), the matrix F G has the 

form F = {Gjji)jji=i^„,^rn, where each block Gjj' G M{n{n + 1)/2,M) is described in Proposition 

3.3. 
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4 A wavelet-based estimator for bivariate OFBM 


In this section, we switch to the bivariate framework (2.12), i.e., n = 2. We draw upon explicit 
expressions for eigenvalues to establish the consistency and asymptotic normality of the estimators 
(1.7) as the wavelet scale grows according to a factor a(z/) —)• oo, as i/ —)• oo. We also show the 
consistency and asymptotic normality, in a sense to be defined, of a sequence of eigenvectors 
associated with the smallest eigenvalue of the sample wavelet variance matrix. 

The proposed estimators make use of the behavior at coarse scales of the sample wavelet 
variance 


1 


K„, 


W{a{u)2^) = y D{a{iy)2^,k)D{a{iy)2^,ky 


a{u)2r 


(4.1) 

fc=i 

Recall that hi < /i 2 under (2.12). In (4.1), {a(i/)}i,gis} is assumed to be a dyadic sequence such 
that 

i' a{v) V 


aiv) < —, 


+ 


0 , V 


oo 


(4.2) 


u a(i/)i+2^i 

(see Remark 4.2 on the choice of a{v) in practice). 

We will make use of some basic relations for bivariate symmetric positive semidefinite matrices. 
Recall that for a matrix 

a b 
b c 

the eigenvalues can be expressed in closed form as 


5 = 


(4.3) 


Ai = 


+ c) - a/A 


< A2 = 


(a + c) + \/A 


A = (a + c)^ — 4(ac — 6 ^). 


(4.4) 


2 - 2 

By Sylvester’s criterion, positive semidefiniteness implies that a + c > 0. As a consequence, if 
det(5) > 0 , 


Ai = ^(a + c)( 1 - 


and 


r 4(ac - 62 )n ^ 1 4(ac - 6 ^) (a+c) 

y (a + c)2 / 2 o + c 4(ac — 62)/(a + ( 

A 2 = 2 + '^) (^ + “ 


4(ac — 6 ^) 
(a + c)2 


(4.5) 


(4.6) 


Let V = (^ 1 ,^ 2 )* G be an eigenvector associated with an eigenvalue A. By further assuming 
6/0, the relation Sv = Av yields 

A — a 


V2 = 


Vl. 


(4.7) 


4.1 The weak limit of eigenvalues 

The next definition describes the proposed estimators for the Hurst eigenvalues /ii < /i 2 - 

Definition 4.1 Let Bh = {B}{it)}t£R be an OFBM under the assumptions (OFBMl)-(4). For 
a dyadic number a(z/), let W{a{h')2y be the associated (symmetric) sample wavelet spectrum at 
scale a{v)2y and let 

Ai(a(i/)2/ < A2(a(i/)2/ (4.8) 

be its eigenvalues. The wavelet estimators at scale a{u)2^ of the eigenvalues hi < 6-2 are defined, 
respectively, as in expression (1.7) with a{v)2^ in place of 2K 
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By analogy to (4.8) and (1.7), we denote the eigenvalues and normalized log-eigenvalues of 
EW{a{i')2^), respectively, by 


Xf {a{u)2^) < X 2 {a{u)2^), hf{a{u)2^) 


logXf{a{u)2^) 

21 og(a(z/)2l) 


< /rf (a(i/) 2 ^) 


logAf(a(;^) 2 J) 
21 og(a(z/)2l) ■ ^ ^ 


When developing asymptotics for the estimators (4.8), in view of the operator self-similarity 
property {P6) we will consider the matrix statistics 

Ba{2^) := p-^Wa{2^){P*)-\ j=ji<...< jm, (4.10) 

where VFa(2'^) = J2k=i j k)D{2^, k)* (c.f. (1.6) and (4.1)). Each Ba{2^) is only a pseudo- 
estimator of 

B{2^) := p-^EW{2^){P*)-\ P G GL(2,M), (4.11) 

because its expression involves the unknown matrix parameter P. It will be convenient to describe 
the matrices (4.10), (4.11) entry-wise as 




(2^)) B{2^) 

\ / li ,22 = 1,2 



^1 ,^2 = 1,2 


(4.12) 


The asymptotic distribution of the matrix statistics (4.10) is given in the following lemma. 

Lemma 4.1 For m G N, let ji < ... < jm be a set of fixed octaves j. Let 11 = ( tt*. ) = 

V ’ / 11,12 = 1,2 

P~^, and let Ba{2^), B{2^) be as in (4.10), (4.11), respectively. Under the assumptions (OFBMl) 
- (4) and (4.2), 


( y^,{vecsBa{2^)-YecsB{2h))) _ 4 iV(0, ^^(ji,..., J,n)), ^ oo. (4.13) 

In (4.13), 

• • • ,im) = di&g^{V)Fdisig^{V)*, (4.14) 

where F is the asymptotic covariance matrix in (3.12), diag^('P) is as defined in (2.1), and 


( T^h 27rii7ri2 tii2 

7rii7r2l 7rii7r22 -f T^12'^21 T^12T^22 

vrli 271217122 7122 


Proof: For any j, a brief calculation shows that vec5(nWa(2-^)n*) = Pvec5lTa(2'^). Likewise, 
vec5(nEVP(2'^ )n*) = V\ecsEW{2^). Therefore, we can recast the left-hand side of (4.13) as 

diag^(P)( y]?“(vec 5 lTa( 2 ^)-vec 5 EW( 2 ^)) 

The weak limit in (4.13) is now a consequence of (4.2) and Theorem 3.1. □ 


The next lemma contains expressions for the wavelet spectrum and its sample counterpart. 
The notation (OiWivdm designates a vector of 2 x 2 matrices. 
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Lemma 4.2 For m G N, let ji < ... < jm be a set of fixed scales j. Let Ba{2^), B{2^) be as in 
(4.10), (4.11), respectively, and let bi.^^ifi2^), bi.^^^ifi2^) be as in (4.12). For a dyadic number a{v), 
under the assumptions (OFBMl) - (f), we can express 


a{v)2^ 




(4.15) 


J —Jl i‘"i3m 


where 


aj,a(F = Pubii{‘^^)a{’^f’^^ + 2 piipi 26 i 2 ( 2 ^)a(i/)'*i +''2 +plfi)22{2^)a{nf^^, 


bj,a{u) = PiiP2\bii{2^)a{vf^^ + {piiP22 + Pi2P2i)bi2{‘2^)a{n)^^"^^‘^ + Pi24'22&22(2^)a(z/)^^^ 

Cj,a{u) = P2]bii{‘2^^)a{yf^^ + 2 p 2 iP 22 &i 2 ( 2 ^)a(i/)'*i +'*2 +plfi)22{2^)a{vf^'^. (4.16) 


Likewise, 


where 


EW{a{v)2^) = ( 

bj,a{u) ^j,a{u) 


aj,a{v) = plibii{2^)a{vf'^^ + 2 piipi 2 bi 2 { 2 ^)a{vf^"^^^ + pl 2 b 22 { 2 ^)a{nf^fi 


(4.17) 


bj,a{u) = PuP2ibii{2^)a{uf’"^ + {puP22 + Pi2P2i)bi2{‘2^)a{n)’"^^’"'^ +Pi2P22b22{2.^)a{vf^^, 

Cj,a{u) = P2ibn{2^)a{i2f^^ + 2p2ip22bi2{2^)a{vf^^'^^ + P22^22(2^)a(jv)2^^ (4.18) 

Proof: The operator scaling properties (P5) and (P6) yield 


W{a{v)2^ 


= [Pdi&g{a{vffia{uf^)Ba{2^)di&g{a{vffia{vf^)P*] 


\hi 


\h2 


hi 


\h2 \ p* 


3 =Jl,..., 3 m ^ ^ 3 =Jl,..., 3 n 

EW{a{v)2^) = Pdiag(a(i/)^^, a(z/)^^) B{2^) (X\a.g{a{v)^^ ,a{v)^'^)P*, 
from which (4.16) and (4.18) follow. □ 


(4.19) 

(4.20) 


The next theorem establishes the consistency and asymptotic normality of the estimators 
described in Definition 4.1. Intuitively, the theorem states that 

hi{a{u)2^) PS hi, h 2 {a{u) 2 ^) ps /12 with convergence rate 2log(a(z/)2-^) 

Theorem 4.1 For m G N, let ji < ... < jm be a set of fixed octaves j. Let hi, / 12 , hf, /if be as 
in (1.7) and (4.9). Under the assumptions (OFBMl) - (4) o,nd (4.2), as v ^ 00 , 


(^) 


(ii) 


(hi{a{v)2^),h2{a{v)2^)) ^ (/ii,/i2), (/if (a(i/)2^),/if (a(i/)2^)) (/ii,/i2) 

for j = jl,.. .,jm; 

f 2log{a{u)2^)^/K~\hi{a{i2)2^) - hf{a{u)2^)] \ 

I 2log{a{u)2^)iyKa,j(h2{a{i2)2^) - h^{a{u)2^)] I 


(4.21) 


/^(O, (jl) • • • ) jm))- 


(4.22) 
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In (4.22) 


(4.23) 


^/il,/i2(jl) • • • 1 jm) — Q^b(JI) • • • jjm)Q j 
where S_b(ji, .. •, jm) is the covariance matrix in (4.14), and Q = (Qjj') 
bloek matrix whose bloeks have dimension 2x3 and satisfy 


^22 (2-^) 


Q.. — detB(2ii) det_B(2-?) 


2bi2(2-^4 (bii{V)b22{2^) 


det i?(22) 

1 

fe22(2'3) 


1 

b22(‘2i) 


and Qj,/ = 0, if j / /. 


Proof: Fix an arbitrary j. For notational simplicity we write B = B{2^) = ) 

V 5^2 — 152 

We also drop the subscripts j, n in the expressions (4.16) and (4.18). Recall that the asymptotic 
distribution for (Ba{2^)] is given by (4.13). 

The statement (i) is a consequence of Lemma A.l and Theorem 3.1. In regard to statement (ii), 
the asymptotics will be written out for just one general term indexed j, but the conclusions apply 
to the whole vector comprising the terms associated with j = ji,... ,jm- Since all meaningful 
limits will boil down to a function of the variables that appear on the right-hand side of (4.13), 
they depend on the omitted, fixed octave j. 

Note that 21og(a(i/)2'^)[/ij(a(i/)2'^) — hf {a{v)2d] = log A* (a(z^)2-^) — logXf'{a{i')2^), i = 1,2. 
Wiht probability going to 1, we can reexpress the latter for i = 1 as 

r, ac — b"^ , ac — 6^ 'I 

\ log - log-^- \ 

t a-|-c a-|-cJ 


+{log(424) 
^ (—4)(ac — 62 )/(a-b c )2 (—4)(ac — 62 )/(a + c)^ / 

The mean value theorem will allow us to show that the second term in the sum (4.24) does not 
contribute to the asymptotic limit. Let 


h{x) = {y/l-x- l)(-a:) ^ / 2 (x) = log/i(x). 

By a second order Taylor expansion of the square root function around 1, 


(4.25) 


fiix) = 


1 -l-(VT^-l) 


+ x^O. 


(4.26) 


Therefore, lim 2 ;_^o — i- Thus, the second term in the sum (4.24) can be reexpressed as 


—4(ac — 6 ^) 
(a -b c )2 


-/.( 


—4(ac — 6 ^) 
(a + c )2 


/ 2 (?afec)| 


4(ac — 6 ^) 4(ac — 6 ^) 


(a + c )2 


(a + c)^ 


where the random variable ^abc hes between ^ and ^ • As a consequence of Lemma 

A.l, it satisfies the limit fate —^ 0; hence, /^(^afec) —^ 4 - By expressions (4.19) and (4.20), 

ac — 1? = a{ uf(h^+h 2 ) Wa(2^), ac-b^ = a(j 2 ) 2 (^i+^ 2 ) det EW{2d. (4.27) 


{ac — b'^) iac — b'^) {ac — b‘^) — {ac — b"^) 
{a + c )2 (a -b c )2 (a + c)^ 


+ (ac -P)| 


(a -b c )2 (a -b c )2 . 
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1 


det Wa{2^) - det EW(2^) 


+ detWa{2^) 


a(jy) 2 (/i 2 -^i) I [0{a{u)^^ ^ 2 ) _|_ +^ 12 )^ 22 ]^ 

1 1 


-[0(a(z/)^i ^ 2 ) + (pf 2 +^* 22 )^ 22 ]^ [Op(a(i/)^i ^ 2 ) + (pf 2 +^’ 22 )^ 22 ]^ 

After premultiplication by the rate factor Kaj, by (4.2) the expression (4.28) becomes 

1 r Op{l) 


}. (4.28) 


a(iy)2(/i2-/Ji) I [0(a(z/)^i“^2) + (p2^ _Pp2^)522]^ 
Op{l) 


I 0 . (4.29) 


+ detiya( 2 ^'), 

■[0(a(i/)'*i-'*2) + (p2^ +p 22 ) 622 ] 2 [Op(a(j^)'^i-'^ 2 ) + (p2^ +^2^) 522 ) 2 ! 

The numerators of the first and second terms appearing in the sum on the left-hand side of (4.29) 
are bounded in probability as a consequence of Theorem 3.1 and of applying the mean value 
theorem to the function / 3 (x) = x^. 

We now look at the first term in the sum (4.24). In view of the relations (4.27) and by applying 
the mean value theorem twice, it can be rewritten as 


log 


detW„(2^) 


-log 


det EWi2^) 


Op(a(i/)'*i-'^2 )(pf 2+^22)^22 "" 0 { a { u )^^ ^") + (P?2+^22)^22 

1 ■ 1 r / (7 ' 

-{det IT, (2^) - detSlT(2^)} - {Op( 






+ 


( 4*12 +7'22)(^22 ~ f^22)|- (4.30) 


The random variable lies between detlTa(2-^) and del EW{2^), whereas the random variable 
^{,22 lies between Op (a (i^) ^^ “^^) +(Pi 2 +P 22 ) ^22 and 0(a(p)^i“^2)_|_(p2^_l_p2^)522. Moreover, these 

random variables display asymptotic behavior —)■ dei EW{2^), ^b 22 (Pu +^ 22 )^ 22 , > 00 , 

where the latter limits stem from Theorem 3.1 and Lemma 4.1, respectively. Therefore, after 
premultiplication by the rate factor Kaj , the expression (4.30) is asymptotically equivalent in 
probability to 


det 5 ,( 2 ^) - det S(2^) ^^^622-622 

V det 5 ( 2 ^) 622 ■ 

Let fi{x, y, z) = xz — y'^. By the mean value theorem, 

det,B,(2^') - detil(2^) = VM^b^^,Ch2,^b22){^^csBa{2^) - yecsB{ 2 ^)), 


(4.31) 


(4.32) 


where ^f^i^bn, ^bi 2 J Cb 22 ) V/ 4 ( 611 , 612 , 622 ) as a consequence of Lemma 4.1. Reintroducing j, 

from (4.29), (4.31) and (4.32), we conclude that 

/;f“{logAi(o(02') - logAf(o(i/)2J)) £ - bn{2>)) 


-2 


bi 2 { 2 ^) (nj\ 5 /oiAA I (bii{ 2 ^)b 22 { 2 ^) A n^fb 22 i 2 ^)-b 22 { 2 ^)\ 

dSsM'^ l'"{ ) “ ” I detB(2^) “ 0 1-M»)-1 ■ 

A similar calculation can be developed for X2{a{i')2^) — Af ( 0 (^) 2 -^). The latter can be recast 


as 


log(a + c) - log(a + c) I I log (^1 


4(ac — 6^) 


1 - 


(d + cy 


- log 1 + w 1 - 


4(ac — 6^) 
(a + c )2 


)}. (4.34) 
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By the mean value theorem based on the function /^{x) = log(l + \/l ~ x), after premultiplication 
by the rate factor y/K^j the second term in the sum (4.34) can be reexpressed as 


/5(ac)v^{ 


4(ac — 6^) 

(a + c)2 


(a + c)2 J 


P 1 

since f'^i^abc) “i by (4.29). As for the first term in the sum (4.34), by the mean value 
theorem we can rewrite it with probability going to 1 as 

\og{Op{a{vf^~^^) + (pi 2 +4*22)^22} - log{ 0 (a(jv)'"i-^ 2 ) + (p ?2 + ^22)^22} 

= ) + bi2 +Pl-2)(b22 - 622)}, 

where ^ 1,22 ( 4*12 +^ 22 )^ 22 - As a consequence, and reintroducing j, 

V^{log(A2(a(02^)) - log(Af(a(02>))} £ ■ 0-35) 

The expressions (4.33) and (4.35) for j = ji < ... < jm imply the weak limit (4.22) with 
asymptotic covariance matrix (4.23). □ 


Remark 4.1 At this point, we can shed more light on the apparent paradox anticipated in the 
Introduction: the asymptotic behavior of each entry of the matrices W{a{v)2^) and EW{a{v)2^) 
(see (4.15) and (4.17)) is generally governed by the power law whereas that of the eigen¬ 

values Ai(a(iA)2'^), Af(a(z^)2'^) is dominated by a(iA)^^b For mathematical simplicity, consider 
only the matrix EW{a{i')2^) under the instance A = P in (1.2), where the columns of P are unit 
vectors. The operator scaling property (P5) of the wavelet transform yields 

EW{a{u)2^) = Pdiag{wia{uf’^\w 2 a{iyf^^)P*, (4.36) 

where EW{2^) = Pdiag{wi,W 2 )P*, for wi = wi{2^) > 0, W 2 = W2{2^) > 0. Then, 

Af (a(j^)2-^) = inf u*EW{a{i2)2^)u < Uq EW{a(u)2^) uq = > (4.37) 

ue5i l|uo|r 

where uq = (P*)“^ei/||(P*)“^ei|| and ei is the first Euclidean vector. In other words, Af (a(z^)2'^) 
cannot go to infinity faster than 

Notwithstanding (4.37), it should be stressed that Xf {a{i')2^) does not generally follow an 
exact power law. This can be seen based on the explicit expression for Xf{a{v)2^) (see (4.4)), i.e., 

Af (a(i/)2^') = ^{EW{a{i2)2^)u + EIT(a(i/)2^')22 - A{aii2)2j)} (4.38) 

where A{a{v)2^) = {EW{a{i')2^)u — ETT(a(z/) 2 -^) 22 )^ + 4 (£'IT(a(iA) 2 '^)i 2 )^. Under (4.36), the 

matrix EIU(a(i/) 2 J) = (EW{a{iy)2^)i,,iA can be expressed entry-wise as in (4.18) with 

V ’ / 21,22 = 1,2 

bii{2^) = wi, 612 ( 2 -^) = 0 and 622 ( 2 -^) = W 2 - It is of interest to note that the quality of the 
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approximation \f{a{i')2^) ~ Ca{v)‘^^^ increases with the difference /i 2 — hi. Indeed, by (4.5) and 
a second order Taylor expansion of the function f{x) = ^/x around 1, 




2 det.BIT(2^) 

+ W2 




det EW{2^) 


1 


-1 1 

. 2 2 -|- W2\'^ 


+ o 


1 


Ol- 


' ^i) y i ^ ^ 

The discussion and conclusions presented in this remark can be easily generalized beyond the 
instance (4.36). 


Remark 4.2 In practice, the choice of a{v) involves a statistical compromise. A large value of 
a{v) with respect to v implies that the estimator variance is relatively large and the estimator 
distribution is not very close to Gaussian, but it also yields a relatively small bias. Computational 
experiments suggest that the ratio vla{y')2^ should be no less than 2^. 


4.2 The weak limit of unit eigenvectors 

For given j, log 2 a{u) £ N, consider the relation (4.7) for the eigenvector entries 

v(a(z/)2-^) = (i;i(a(i/)2'^),U2(a(i/)2'^))* G (4.40) 


associated with the smallest eigenvalue Ai ;= Ai(a(j^)2-^) of the (symmetric) sample wavelet 
variance matrix S := W{a{i^)2^). Further consider their wavelet variance counterparts 
(ui(a(z/)2'^),U2(a(z/)2'^))*, Xf := Xf{a{iy)2^) and S := EW{a{iy)2^). The ratios {v2/vi){a{i')2^), 
{v 2 /vi){a{u) 2 ^) represent the tangents of the angles that determine the associated eigenspaces. 
In the next definition, we propose to use {v2/vi){a{i')2^) as an estimator of the tangent of the 
angle determined by the entries of P associated with Ai when P £ 0(2), i.e., tan(0) = —^ 12 /^ 22 - 
For notational simplicity, we will just write 

(4.41, 

P22 

Definition 4.2 Let j, log 2 a(i^) £ N. Under the assumptions of Definition 4.1, we define the 
estimator of 0 (see (4.41)) at scale a(v)2^ and its wavelet spectrum counterpart as 


6(a{v)2^) 



(a{v)2^) 


Xi{a(v)2^) -aj,a(u) 

bj,a{u) 


6(a{i')2^) 



{aiu)2^) 


Xf{a(v)2^) - 

bj,a{u) 


(4.42) 


where the definition of each individual term on the right-hand of the expressions in (4.42) is given 
by (4.8), (4.9), (4.15) and (4.17). 


The consistency and asymptotic normality of the estimator 6(a{v)2^) in (4.42) are precisely 
stated and shown in Theorem 4.2 below. The limits themselves do not depend on the orthogonality 
of P. Note that the parametric assumption (4.43) below rules out diagonal wavelet spectra, the 
latter being associated with entry-wise scaling instances of OFBM (as in (1.4)). For further 
comments on the assumptions of Theorem 4.2, see Remark 4.4. 
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Theorem 4.2 Let j G N, and let 9{a{v)2^), 0{a{v)2^) he as in (4.42). Suppose the assumptions 
(OFBMl)-(4) and (4.2) hold, as well as the parametric assumption P 22 / 0 . 


(i) If 

either pi 2 / 0 or pi 2 = 0 and 612 / 0, (4.43) 

then 

9{a{h')2^) 9 = , 9{a{i')2^) ^ 9 = , u ^ oo; (4.44) 

P22 P22 

(ii) if 

P12 7 ^ 0 and 612 / 0, (4.45) 

then, 

a{vp-^^^/K~{0{a{v)2^)-9{a{u)2^)]AN{Q,al), ^ 00 . (4.46) 

In (4.46), cj| = TZ*Tib{‘^^)'R,, where Sb(2'^) is the 3x3 block, associated with j, on the main 
diagonal of (see (4.13)^, and TZ* = -1,’^y 


Proof: As in the proof of Theorem 4.1, we will drop the subscripts j and 12 whenever convenient. 
To show (i) for 9{a{v)2^), first assume that pi 2 / 0. The property (P7) ensures that 622 / 0, 
whence h = hj^a{v) / 0 for large v in expression (4.18). The same applies to a = aj a[u)- Therefore, 
ui / 0 in (4.7), and 


V2 

Vl 


a/Af 


- 1 


Pl2 
P22 ’ 


V —)• 00 , 


(4.47) 


by Lemma A.l. 

Now assume that pi 2 = 0. Then, pn / 0, and a = aj^a{u) / 0 for large 12 . Since 612 / 0 by 
assumption, then b = hj^a(v) / 0. Again, the limit (4.47) follows. 

In regard to 9{a{v)2^), in either case in (4.43) the relations above can be simply rewritten 

for S = W{a{i2)2^) with eigenvector v. Because of (4.13), Ba{2^) B{2^) and thus the 

expression V 2 IV 1 is well-defined with probability increasing to 1, as —t- 00 . Again by Lemma 

a.i,u2M4-^. 


For the ensuing developments, recall that 622 ^ i® well-defined. So, to show (ii), reexpress 


./K~ 




= n/^{^ - 2}' G®) 

The first term on the right-hand side of (4.48) can be further developed into 

In view of Lemmas A.l and 4.1, 

Af , 1 _ detW(2^') ^h-b^ P ^ {b22-h22} d N{0,a\b22{2l))) 

b a(jv) 2 (/* 2 -/*l)p^ 2 P 22 &i 2 ( 7 ’l 2 +^' 22 )’ ^ 6 -I ^ 622 & 22 ( 2 ^') 
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where (T^( 622 ( 2 -^)) comes from the matrix Sb(ji, • • • ,jm) in (4.13). Consequently, in regard to the 
second term in the sum (4.49), 




(4.50) 


We now look at the first term in the sum (4.49). Let f = 4^j^ and r = 4^^^. Then, 


c—6^ 


V^{Ai-Af} = v^{ 2 ( 

'ac — b'^ 


ac 


-P 


a + c 


- 2 


ac 


- 6%'I /\/l -r- 1 


a + c 


—r 


+2 


a + c 


)v^{(- 


\/l — r — 1\ /\/l — r — 1 


—r 


—r 




(4.51) 


By the same reasoning leading to (4.29), 


VTZi{2{ 

+ det EW(2^) 


ac — b‘^\ ^/ac — b"^ 
a + cJ \ a + c 


)}=o( 0 “'{ 

Op(l) 


Op(l) 


Op(a(p)^i-^ 2 ) + (pfs +^22)^22 

}. (4.52) 


[Op(a(zy)'^i-^2) + (p2^ +p 22 ) 622 ][ 0 (a(z^)'^i-'^ 2 ) + {pj^ +pi2)M 

Moreover, for /i(x) as in (4.25), with probability going to 1 the mean value theorem gives 
\/^aj{hP) ~ /i(^)} = fi{P)\/^a,j{r — r), where the random variable lies between r and f. 
Thus, by (A.3), (A.4), (4.26) and (4.29), 


fUc \ K ^ ^ P Op{l) 


(4.53) 


Therefore, by (4.52), (4.53), (A.3) and (A.4), the first term in (4.49) is asymptotic equivalent in 
probability to 


a(p) 


h2 — h\ 


1 


[Op(a(p)^i ^2)+^^^2^22^22] 


Tj){a(p)->Op(l) + -Tr«^}4o. (4.54) 


As a consequence, in regard to the first term on the right-hand side of (4.48), 


We now turn to the second term on the right-hand side of (4.48). Note that 


T_T|4o. 

b b ) 


a a a{b — b) — b{a — a) p 2 h 2 ^^b — b) — bia — a) 

b b hb ^ ^ Pl2P22b22p^) 


(4.55) 


(4.56) 


Note that the ratio on the right-hand side of (4.56) is well-defined, since pi 2 ,P 22 / 0, by assump¬ 
tion. Recall the expressions (4.16) and (4.18) and denote a{u)~‘^^^b, and 

a{iy)~‘^^^b, respectively, by 

0 (a(i.) 2 (^i-'^ 2 )) + aua{i2p^-^^ + 022 , 0 (a(i/) 2 (^i-^ 2 )) + + ^ 22 , 

Op(a(p) 2 ('^i-'^ 2 )) + ai 2 a{j 2 p^-’^^ + S 22 and Op(a(p) 2 (^i-^ 2 )) + + ^ 22 - (4.57) 
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For the sake of mathematical clarity, we will henceforth omit the terms of order in 

the expressions (4.57), since the ensuing argument can be easily adapted to account for them. So, 
based on (4.57), we can rewrite the numerator in (4.56) as 

+ a22){(/3i2 — + {1322 — /322)} 

— (/3i2a(z/)^^“^^ + /322){(Si2 — ai2)a{v)^^~^'^ + (S 22 — 022 )} 

= {{ai2a{v)^'^~^'^ + a22)(/3i2 — /?12) — {I3i2a{i')^^~^'^ + /322)(Si2 — ai2)} 

+a{v)^^~^^ {ai2{(322 — 1322) — (3i2{a22 — 0122 )}, (4.58) 

where the equality is a consequence of the fact that a 22(/322 — 1322) — l322{ot22 — 022 ) = 0. When 
premultiplied by a(z/)^ 2 -^i ^Kaj, the hrst term in the sum (4.58) is asymptotically equivalent in 
probability to 

^/K^{a22{l3i2 — /3 i2) — I322{oci2 — 012 )} —^ —7'i 2&22(2'^) det(P) A^(0, cj^(6i2(2'^))), (4.59) 

where cj^(6i2(2'^)) is taken from the matrix Sb(ji, • • • , jm ) in (4.13). In turn, when premultiplied 
by a{v)^‘^~^^ the second term in the sum (4.58) is asymptotically equivalent in probability 

to 

^/i^j{al 2{^22 - (322) - I3i2{a22 - 022 )} 4 p? 2 & 12 ( 2 ^) det(P) iV(0, ( 622 ( 2 '))). (4.60) 

In view of (4.55), and assuming P 22 > 0, the relation (4.46) is obtained by dividing the sum of 
the weak limits (4.59) and (4.60) by the denominator on the right-hand side of (4.56). □ 


Remark 4.3 From a different but mathematically equivalent perspective, the statement ( 4 . 44 ) 
implies that W{a{i2)2^) and EW{a{i') 2 ^) have sequences of unit eigenvectors associated with 
Xi{a{i') 2 ^) and Xf {a{i') 2 ^), respectively, which converge (in probability and deterministically, 
respectively) to the limiting unit vector 




IP22I -Pi2 sign(p22) \ 

+ Ph ’ Vpi2+P22 ^ 


( 4 . 61 ) 


If P G 0 ( 2 ), then eigenvectors of H in two orthogonal directions can be consistently estimated 
by the eigenvectors of W{a{v) 2 ^). In view of ( 4 . 61 ), there is also a unit eigenvector associated 
with X2{a{u)2^) that converges to (p2i, —Pii)*, since the latter is orthogonal to {P22, —P12)*■ It 
is easy to see that the latter vectors generate the same (eigen)spaces as (^12,7*22)*, (pii,7'2i)*, 
respectively. 


Remark 4.4 In regard to the assumptions of Theorem 4 . 2 , when P22 = 0 , it is clear that ( 4 . 42 ) 
cannot be consistent. By a similar proof to that of Theorem 4 . 2 , asymptotic normality (with 
different variances) for {6{a{i2)2^) — 9{a{i2)2^)} can also be established when pi2 = 0 or 612 = 0. 
In the former case, for instance, the convergence rate is a{v)^'^~^^ or respectively, 

depending on whether 612 / 0 or both 612 = 0 and p2i / 0. 

Remark 4.5 In Theorem 4 . 2 , the convergence rate is the non-standard a{i')^^~^^ Kaj ~ 
a(z/)^2-^i-i/2^jyy'2i^ which depends on the parameters to be estimated hi, 6.2. In practice, since 
a(z/) is much slower than 1/ (see condition ( 4 . 2 )), its effect may not be noticeable (see Section 5 
below). 
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Remark 4.6 Substituting X2{a{i')2^) for Xi{a{iy)2^) in (4.42) does not provide information on 
pii or P 21 - Under mild assumptions on the parameters, Lemma A.l yields 

g / A 2 (a(i^) 2 -^') P22 

b\ a ) pi 2 


5 Simulation studies 

In this section, we lay out and discuss a broad computational study of the performance of the 
proposed estimators. OFBM was simulated by means of the Hermir Toolbox, devised in Helgason 
et al. (20116, 2011o) and available at www.herniir.org. When describing the results, we drop 
the asymptotic scaling factor a(p) used in Theorems 4.1 and 4.2 and only speak of shifting scales 
2^ e N. 

The chosen sample path size was N = 2^®. The purpose of picking a large N was to provide a 
compelling illustration of the estimators’ ability to capture both Hurst eigenvalues. In biomedical 
applications, typical recordings can be much shorter (of the order N = 2^^; see, for instance, 
Ivanov (1999)). Nevertheless, sample paths of size N = 2^® are, indeed, encountered in Internet 
traffic analysis (Abry et al. (2002)) and hydrodynamic turbulence (Frisch (1995)). 

In the wavelet analysis, we used least asymmetric orthogonal Daubechies wavelets with = 2 
(Daubechies (1992), chapters 6-8). It was verified that similar conclusions can be obtained when 
using > 2 or some other wavelet with a large enough number of vanishing moments. 

Remark 5.1 In practice, a continuous time OFBM sample path is not available and thus the 
theoretical wavelet coefficient D{2^, k) cannot be computed. Instead, one approximates the latter 
by means of the classical recursive (or pyramidal) discrete filter bank algorithm (Mallat (1999), 
chapter 7). The algorithm’s main input is a discrete time sequence, which, following Veitch et al. 
(2003), can be a discrete OFBM sample {BH{k)}k£T, T C Z. In Section C, we lay out in detail 
the mathematical framework for estimation based on discretized wavelet coefficients. In the main 
result of the section. Theorem C.l, we establish that, under mild conditions, the estimated Hurst 
eigenvalues and eigenvector stemming from the pyramidal algorithm also satisfy the weak limits 
(4.22) and (4.46), respectively. 


5.1 Entry-wise vs eigenvalue-based estimation 


A matrix IU(2'^) was computed from a single sample path (of size N = 

0.98 0.57 Y ^ diag(0.25,0.85). 


with matrix parameters P = 


0.20 0.82 


2^®) of a synthetic OFBM 
In the spirit of the entry- 


wise approach, the wavelet-based estimation of the Hurst eigenvalues relies on performing a linear 
regression on a log 2 W(2^) vs j = log 2 2^ diagram, motivated by the log-transformed scalar version 
of EW{2^) as in (H5), i.e., EW{2^) = C{H)2^'^^ for some C{H) > 0. This is shown for each 
auto- and cross-wavelet components of the bivariate OFBM in Figure 1. The top row displays, 
in order, plots for log 2 IU(2-^)i^i vs j = log 2 2^, log 2 IU( 2 '^)i ^2 vs j = log 2 2^ and log 2 IU( 2 '^) 2,2 vs 
j = log 2 2P The asymptotic behavior j2/i2 is superimposed on each of these plots. In view of the 
expression (4.16), it is unsurprising that at coarse scales all auto- and cross-components end up 
driven by the dominant Hurst eigenvalue /i 2 . In other words, the conspicuous prevalence of the 
latter precludes the estimation of the Hurst eigenvalue hi. By contrast. Figure 1, bottom row, 
displays, in order, plots for log 2 Ai(2'^) vs j = log 2 2^ and log 2 X 2 { 2 ^) vs j = log 2 2^, as well as the 
superimposed asymptotic behaviors j2hi and j2/i2, respectively. The eigenvalue-based procedure 
leads to two scaling laws individually driven by each of the Hurst eigenvalues hi and /i 2 . 
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5.2 Estimation performance and asymptotic normality 

We also numerically synthesized 10, 000 bivariate OFBM sample paths to study the finite-sample 
effectiveness of the normal approximation described in Theorem 4.1. For each path, the estimates 
W{2^), Ai(2-^), A 2 ( 2 -^) and V2{‘2^)/vi{‘2^) were computed. Averaging over realizations yields Monte 
Carlo estimates .EAi( 2 '^), E\2{2^) and £^pi 2 ( 2 '^)/p 22 ( 2 -^) of the ensemble averages E\i{2^)^ E\2{2^) 
and —E9{2^), together with estimates of the variances VarAi(2'^), VarA 2 ( 2 -^), Var0(2'^). 

In Figures 2, 3 and 4, the simulated OFBM has parameter Jh = diag(0.25,0.85) and 
path size N = 2^® for the three cases, but different mixing matrices of the general form 

p — ( V In this parametrization, (3 = —9 (see (4.41)). The simu- 

\7/vl+7^ 1 /Vl + W y 

lated instances are representative of different parametric settings. The choice (3 = 0.7, 7 = 0.2 
(Figure 2) illustrates the situation of a general mixing matrix P. The choice (3 = —7 and 
P /\/l + /3^ = sinyr/G (Figure 3) represents the case P £ 0(2), whereas 7 = 0 and /3 = 0.2 (Figure 
4) portrays the case often referred to as fractional cointegration. A comparison of the estimates 
(black solid lines with ‘ 0 ’) with the theoretical values (red dashed lines) in Figures 2, 3, and 4 
reveals the robustness of the proposed estimators (top row). The qq-plots (bottom row) also show 
that no deviation from a AA(0,1) distribution can be observed within ±2 standard deviations for 
log 2 Ai( 2 '^) and log 2 Ai( 2 '^), and within ±2 standard deviations for pi 2 fp 22 - 

The simulations indicate that, beyond asymptotics. Theorem 4.1 for hi{2^) and h2{2^), and 
Theorem 4.2 for 9{2^) = —/3(2-^) provide effective normal approximations to the finite-sample esti¬ 
mator distributions. Whether or not P £ 0(2) does not impact the performance of log 2 Ai(2'^)/2j, 
log 2 A 2 ( 2 '^)/ 2 j and (3(2^). However, assuming P £ 0(2) additionally allows for the full estimation 
of P, and thus of PI, as discussed in Remark 4.3 and illustrated in Figure 5. 


5.3 Beyond the bivariate setting 


It is natural to ask whether in higher dimension the eigenvalues of the sample wavelet variance 
are good estimators of the Hurst eigenvalues. Figure 6 shows eigenvalue-based estimation at work 
for n = 4, with log 2 Xp{2^)/2j, for p = 1,2, 3,4, averaged over 2,000 realizations of an OFBM with 
parameters 


P = 


( 0.90 
0.43 
0 

V 0 


-0.22 

-0.30 

-0.22 \ 


0.45 

0.63 

0.46 


-0.85 

0.40 

0.30 


0 

-0.59 

0.81 y 


= 2^®. 

The computational : 


and Jh = diag(0.20,0.40,0.70,0.90), N 
smallest and largest eigenvalues of W{2^) are still good estimators of the smallest and largest 
entries of Jh- Furthermore, there is evidence that the method produces reasonable estimates of 
all the intermediate-valued entries of Jh based on the corresponding intermediate eigenvalues of 
W{2^). In regard to eigenvector estimation under the assumption P £ 0(4), unreported numerical 
studies suggest that the sample wavelet variance eigenvectors are also good estimators for P. 


6 Perspectives and open issues 

OFBM constitutes the natural multivariate extension of the univariate FBM, allowing a different 
Hurst eigenvalue in each coordinate, in an arbitrary coordinate system. When the mixing matrix 
P is non-diagonal, the problem of estimating H becomes distinctively multivariate and is not 
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Figure 2; Estimation performance and asymptotic normality: OFBM with 7 = 0.2 and 

(5 = 0.7. Top row: black solid lines with ‘o’ show ± \J Var 'd/n for -d = log 2 Ai(2'^ )/2j, 

log 2 A 2 ( 2 -^)/ 2 j and pi 2 ( 2 '^)/p 22 ( 2 -^) (target parameters: hi (left plots), /i 2 (center plots), P 12 /P 22 
(right plots), respectively), red dashed lines corresponding to theoretical values. Bottom row, the 
corresponding qq-plots (against AA(0,1) distributions) for j = 10. 


easily amenable to approaches inspired in the univariate context. In this work, we propose a 
change of perspective from univariate-like, entry-wise scaling relations to the eigenstructure of 
the sample wavelet variance W{2^) across scales. In the bivariate setting, this methodology 
is mathematically shown to yield consistent and asymptotically normal estimates of the Hurst 
eigenvalues as well as of the eigenspace angle parameter —pi 2 /p 22 - Consequently, the matrix 
Hurst parameter H can be fully estimated when P G 0(2). Large sample size simulation was 
used to illustrate and shed further light on the weak limits obtained. The research contained 
in this paper has lead to five open issues, currently under investigation: (i) the quantitative 
assessment of the performance of the estimators as a function of sample size; how much data 
is demanded by the difficult problem of estimating operator-scaling systems, especially in high 
dimension?; (ii) is there an advantage to nsing multiple scales in a regression, as in univariate 
wavelet-based estimation?; {iii) mathematical extensions to OFBM in higher dimension; (iv) the 
estimation of non-orthogonal coordinate systems; (n) applications in real data. In the near future, 
a Matlab toolbox for the estimators proposed in this paper will be made publicly available. 


A Additional results 


Proof of Proposition 3.1 We first show (PI). Since the covariance function EBH{s)BH{t)* 
is continuous, by Cramer and Leadbetter (1967), p. 86 , it suffices to show that 


\\EBH{s)BH{t)*\\ii \'ip{s)\\'ijj{t)\ dsdt < 00 . 


(A.l) 
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Figure 3: Estimation performance and asymptotic normality: OFBM with 7 = —j5 

and /3/v^l + = sinvr/h. Top row: black solid lines with ‘o’ show E-d ± \J Var 'd/n for 

"d = log 2 Ai( 2 -^)/ 2 j, log 2 A 2 ( 2 '^)/ 2 j and pi 2 ( 2 -^)/p 22 ( 2 '^) (target parameters: hi (left plots), /i 2 
(center plots), P 12 IP 22 (right plots), respectively), red dashed lines correspond to theoretical val¬ 
ues. Bottom row, corresponding qq-plots (against AA(0,1) distributions) for j = 10. 


In fact. 


\\EBH{s)BHm\p = E E \EBH{sh,BH{t)i,\ < ( E \lEB h{s)1) ( E • 

2 l = l 22 = 1 2 l=l 22 = 1 

(A.2) 

However, for t E M, \\EBH{t)BH{t)*\\ioo = H^oo < (7|i|2max5Reig(/f)_ Therefore, (A.2) is 

bounded by (7|i|™a'x5Reig(H)|g|max5Reig(H)_ gy conditions (2.7) and (2.8), (A.l) holds. The fact that 
E\\BH{t)\\ii < 00 and the assumption (2.7) yield 

E [ \m\\\BH{2H + 2^k)ydt<C [ \m\ EWBnit + k)\\iidt 

< CE\\Bh{1)\\ [ |V'('S - ^)|||s-^||pds < 00 , 

Jr 

where we used the change-of-variables s = t + k. By Fubini, this yields (PI), as claimed. 

In view of (PI), the properties (P2), (P3) and (P5) and can be established by arguments 
similar to those for the univariate case (see, for instance, the argument in Delbeke and Abry 
(2000), Theorem 3 for the former two, and Bardet (2000), Proposition II. 1, for the latter). Note 
that = C^p{2^x) for some constant C > 0. Therefore, the property (P6) is a 

consequence of (P3). The property (P4) is a consequence of the harmonizable representation 

(1.2) , and also of the conditions (2.6), (2.7) and (2.8); the latter ensure that the integrand in 

(3.2) is well-defined in M. We now show (P7). Since 'ijj then 'tjj(—2^x) = ijj{2^x). Thus, by a 
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Figure 4; Estimation performance and asymptotic normality: OFBM with 7 = 0 and (5 = 

0.2. Top row: black ‘o’ in solid lines show E'd ± for 1 ? = log 2 \i{2^)/2j^ log 2 \2{2^)l2j 

and Pi 2 ( 2 -^)/ 4 ' 22 ( 2 -^) (target parameters: hi (left plots), /i 2 (center plots), P 12 /P 22 (right plots), 
respectively), red dashed lines correspond to theoretical values. Bottom row, corresponding qq- 
plots (against AA(0,1) distributions) for j = 10. 

change-of-variables y = —x over the integration domain x < 0 we can rewrite (3.2) as 

EW{2^) = C2^ [ x-^2?R:{AA*)x-^*&^^fi^dx. 

Jo ^ 

Since supp ij^ix) has positive Lebesgue measure, then (2.10) yields v*EW{2^)v > 0, u G C"\{0}. □ 


Lemma A.l Fix j G N. Under (OFBMl)-(4) and (4.2), let EW{a{u)2^) and W{a{i')2^) be as 
in Lemma 4-2. Then, the following limits hold, as n ^ 00 : 


^ac-b^) 4detEW{2^) 2/.1 

« + c ^ {Pj2+Pl2)b22{2^f^’'^ 


4(ac-62) 4detElF(2^) 1 

(a + c)2 ~ {{pj,+pl,)b22{2^Wa{nyi'^^-'^U > 


d{ac-P) p 4detlT(2J) ^/A 2 fei 

® + ^ (Pi 2 +^ 22 )^ 22 ( 2 ^') 


4 (aE-P) P 4detlF(2^) 1 

(a + c )2 ^ {{pl2+pl2)b22{2^WP’^?^^^-^^'>' 


4(ac — b‘^)/{a + c)^ 
det EW{2^) 


1 

2 ’ 


(1 . /i 4(ac-b2)\ 

p 1 

4{ac — h"^) / {a + cY 2’ 


{Pi2+P22)b22{2Y ^ 

^2 ~ {Pi2 + P22)b22{2Ya{i'f^Y 


I P 
Ai ~ 


W 2 +Pi 2 )J> 22 (») 

A2 ~ {pu +P22)^22(2^)a(l/)^^T 


(A.5) 

(A. 6 ) 

(A.7) 
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Figure 5: Estimation performance when P G 0(2): OFBM with 7 = —f3 and (5/^/1 + [3'^ = 
sinvr/G. Estimates k,l G {1,2}^ (black solid lines with ‘ 0 ’) of the entries of P, pk^i (red 

dashed line), from the eigenvectors of the W{2^). When P G 0(2), both P and Jh can be 
estimated, and the matrix exponent H is fully estimated. 

Proof: In regard to Ai and Xf, consider the relation (4.5). By Lemma 4.2, 

4(ac — }p‘) 
a + c 

_ _ 4a(i/)^(^i+^^) det EW{2^) _ 

(Pll + + 2(P11P12 + P2lP22)bl2a{l')^^+^'^ + {p\2 + P22)^22a(l^)2^2 ’ 

where the determinant is non-trivial due to the property (P7). Again from Lemma 4.2 and by 
applying Theorem 3.1, an analogous expression holds for 4(ac —6^)/(a + c). Then, the expressions 
(A.3), (A.4), (A.5), (A.6) follow. 

An analogous reasoning applied to A 2 and Af leads to (A.7), since the relation (4.6), Lemma 
4.2 and Theorem 3.1 show that Al' ~ a + c, A 2 ~ a + c. □ 


B Asymptotic normality at fixed scales: proofs 

This section contains the remaining proofs of the statements in Section 3. 


Proof of Proposition 3.2: In view of (2.7), we can assume without loss of generality that 

supp(V^) = [0, iL], K>1. (B-1) 

By (2.6), (2.11), (2.13), (3.1) and (B.l), the wavelet covariance can be reexpressed as 

ED{2\k)D{2^',k')* = [ [ i>it)^lj{t')EBHi2H + 2^k)BH{2^'t'+ 2^'kydtdt' 


^p{t)^p{t')\2H - 2^'t' + 2^k - 2^'k'fE\2H - t' + 2^k - 2^'k'f”dtdt' (B.2) 






1 + 


2H - 2^'t' 


2ik - 2^'k' 


H 


1 + 


2H - 2^'t' 


2^k - 2^'k' 


H* 


dtdt 


Let 


f{x) = = ifi^,i 2 {x)) 

\ / 21,22 = 1,.••,71 


'^\2^k-2^' k'f\ 
(B.3) 

(B.4) 
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Figure 6 : Estimation performance and asymptotic normality: OFBM in dimension 4 . 

Top row: black solid lines with ‘o’ show E-d ± for i? = log 2 Ai(2-^)/2j, log 2 A2(2-^)/2j, 

log 2 A 3 ( 2 -^)/ 2 j and log 2 A 4 ( 2 '^)/ 2 j (target parameters: /ii, /i 2 , hs, h^, respectively, with the plots 
in the same order), red dashed lines correspond to theoretical values. Bottom row, corresponding 
qq-plots (against AA(0,1) distributions) for j = 10. 


Within the radius dehned by (3.4), the integrand in (B.3) can be reexpressed analytically as 


1 + X7 


2H-2U'\H / 2H-2^'t'\H 


2ik-2^'k') 


S 1 + — 


21k- 21 k' 


r\ \2ik- 


r=0 


2-^ h/ / 


Thus, again (2.6) yields 

pK rK 

Jo 


JO 


^ 2it - 2i't' \ H / 211 - 2i't' . 

1 H- r-;-r— ) S ( 1 H- r-- Tj— ) dtdt 

-2ik'J V 21k-21k') 

/2if — 2i'f'\r \ 

I {:2ik-2i-k’) .„■ 


CX3 . 

h,i2 


E 

r=2N^ 


21k 


r! 


'0 JO 


We now look at each term in the summation (B.5). By induction, up to the matrix constant 
(/{2N^ + /)!'j the term associated with the index r = 2N^ + 1, I £ N, can be 


\lh,i2 

written as 


h,*2=1, 


{2ik - 2i'k'f^'i I f 




E 


v=N, 


V' 


-V rK rK 


2N^ +1\ (2^y(-2J'')(^^^+^) 

V ) {2ik — 2i'k'Y Jq Jq 

t=u-N^ {21^1')^'!’ 

{2ik - 2i'k'f^^ 


dtdt'^ 
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{2jk - 2i'k'y 




dtdt'Y 


In the expression (B.6), for any I G N the largest binomial coefficient is 


(B.6) 


f 2N^ + l \ l_ 

\\i/2\+nJ^ yi’ 


(B.7) 


where the asymptotic equivalence (as I —)• oo) follows from Stirling’s formula. Moreover, 


rK t-K 




< C 


i ■ 


(B.8) 


Recall that we are taking the wavelet parameters in the range (3.4). Then, by (B.7), (B.8) and 
Lemma B.l below, the expression (B.5) is bounded from above by 


||/(2Ar^)(l)|| /2ArA 

{2N^)\ \N^J 



\i;{t)i;{t')\t^^{tY^dtdt' 


^ ||/(2A^^+0(1)|| max{2^2^y d,^+ 

^ (2iV^ + 0! 2ik - 2Tfc' I + Jo 





1=1 ^ ^ 




Therefore, in view of (3.4), the series (B.5) is, indeed, summable, and by (B.3) and (B.6) the 
expression (3.5) holds. 

In regard to (3.6), note that for c > (5“^Rrmax{2'^ , 2-^ } we can construct the bound for the 
matrix exponential 


llc^ll < C ||P||||P“^||||y"||p < C" c”'"^h6eig(H) 5RW |log"’(c) 
Now set c = \2^k — 2^ k'\ under the restriction (3.4). □ 


Lemma B.l Let f be as in (B.4). Then, 


||/«(l)||/(r-l)! = 0(I), reN. 


(B.9) 


Proof: To prove (B.9), we will first establish the formula 


(r-l)-fc 


k=0 ^ ^ ji=0 


H*-kI 


k-l 

n( 77 * 

32=0 


hi)}, 


(B.IO) 


where, for notational convenience, we set nj=o(-^ ~ jl) ^ I- r = 1, the formula holds. So, 
by induction, assume that it also holds for r G N. Then, 


/(^+i)(x) = 


r 

= E 

fc =0 


;){ 


(r-l)-fc 

n i^-hi 

31=0 


H-(r- k)I)x^hr-k)l-l^y. 


H*-kI 
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k-1 

- A:/)}| ll {H*-j2l)] 

32=0 

V / j=o 

+E [([) + (j! i)]{ n(if-*/)} 

A:=l ^ ^ ^ ^ ji=0 j2=0 

+ - in} 

^ ^ 3=0 

/y* I \ '^ ^ ^ 1 

fc=o ^ ii=o i2=o 

This establishes (B.IO). 

Recall the notation (2.9) for the eigenvalues of H, and let /imin := minifteig(ff) G (0,1). For 
any J G N U {0}, the equivalence of matrix norms yields the bound 

J J J n j 

jl[{H-jI)\\<Cll\\JH-jiy=Cll{Y,\j-hk\+c}<C']J\j-h^i^\, (B.ll) 

j=0 j=0 j=0 k=l j=0 

where c > 0 accounts for the number of off-diagonal Is in the Jordan form Jh- Now set x = 1 in 
the expression (B.IO). By (B.ll), 

||/M(i)|| = |(-i)n||^Q{ n 

k=0 ^ ^ ji=0 j2=0 

n IR “ |j2 - /iminlj = C'|/^J.^(1)|, 

k=0 ^ ' ji=0 j2=0 

where /vin(^) denotes the function / when n = 1 (scalar), the Hurst parameter is /imin and 
S = 1. Since = ‘^h.amlVjZlij - 2/;-min), then (B.9) holds. □ 


For the next results, we assume that aj, aji G N. 

Let {4>k}k£Z be a sequence of real numbers. We are interested in calculating the limit 


^ Kj Kj' 

lim — \ ^ \ ^ 4*aik—a,ik'- 
k=l k'=l 

(B.12) 

TZ := {ajfZ — Uj/N} 

(B.13) 


be the set of indices of (j).. The next lemma provides the growth rate of the number of terms in 
the summation (B.12). Before we state and show it, note that to every r £ TZ we can associate 
the affine level curve 

k\k) = ^k - keR, (B.14) 


Qj'j^ (X')f 


which contains all the pairs (/c, k') G satisfying ajk — aj/k' = r in the chosen range. 
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Lemma B.2 Let aj,aj>,i' G N and r G TZ (see (B.13)J. Let k'{-) he the function defined in 
(B.14). Then, the number of pairs {k,k') G 1?, 1 <k < aj>n, 1 < k' < aju, in the affine 
level curve (B.14) associated with r satisfies 

lim ^ = gcd(a,', a,'). (B.15) 

u^oo 1/ 

Proof: We first show that ajk — aj/k' = r has a solution (/c*, k'{kfij) G 1?, and that the set A 
of such solutions has the form 

A = { (/c**, /c'(A:**)) G : A:** = /c* H- -r^ -rZ|. (B.16) 

t gcd(aj,aj/j ) 

By Theorem 1.8 in Jones and Jones (1998), p.lO, we can conveniently reexpress the set TZ in 
(B.13) as 

TZ = gcd(aj, aj')Z. (B.17) 

In other words, for any fixed w G Z, there is a solution pair (A:*, k'{kfij) G for the equation 


Ujk 


Oj! k' 


gcd(aj, aji)w. 


(B.18) 


For notational simplicity, rewrite the solution pair as (A:*,A:(), and consider another solution 
{k^*,k(fi} G Z?. By plugging the two solution pairs into (B.18) we obtain 


A:** ~ A:* 


aj/gcd(aj,aj/) 

/ 1 / X I 

ay j gcd(aj, ay) 


(B.19) 


By premultiplying (B.19) by —1 if necessary, we can assume that A:** — A:*,A:(^ — A:( G N. Since 
Ojj gcd(aj, Oj/), ay j gcd(aj, a^/) are coprime, then A:** —A:* = gcd(a. a-,) ^ some z G N. Therefore, 

(A:**, A:(^) has the form described on the right-hand side of (B.16). Conversely, for a given solution 
(A:*, k'{kfi)), every pair (A:**, A:'(A:**)) of the form expressed on the right-hand side of (B.16) is also 
a solution (in Z^) for (B.18). This establishes (B.16). 

To show (B.15), for any sufficiently large v, let ko G {1,..., aj>n} be the smallest number 
such that {ko,k'{ko)) G solves (B.18). Let M. 3 x = gcd{aj,aj'){i' — ko/ajfi. Then, by (B.16), 
([xj, A:'([a:J)) is the rightmost solution for (B.18) within the first-entry range k = 1,... ,ajin. 
Moreover, —)■ gcd(aj,aj/), n -3 oo. This establishes (B.15). □ 


Lemma B.3 Let {fi.} G M 6e o sequence such that Y((T=-oo \'i^zgcA(aj,a-,)\ < oo. Then, 

ayu ajU oo 

fiajk-ayk' ^ gcd{aj,ay) ^ (l>zgcd{aj,ay), V ^ OO. (B.20) 

k=l k'=l z=—oo 

Proof: By expression (B.17) in the proof of Lemma B.2, X]a:=i = 

where Bjjfiv) is the range for r such that 1 < A: < ajiu and 
1 < A:' < Ojn. Moreover, by (B.17), U^izNBjyfin) = gcd(aj,aji)Z. Then, expression (B.20) is 
a consequence of dominated convergence and (B.15). □ 
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Proof of Proposition 3.3: The statement (ii) is a direct consequence of (z), so we only prove 
the latter. It suffices to consider the subsequence p = 2^2^ p*, p* —)■ oo. Then, Kj = 2^ u* , 
Kji = and ^JlZjKJ^KJ,^ = . From the expression (B.2), 

CoY{D{2^,k),D{2^',k')) =: 

Let '^ 2 ik- 2 i'k' ~ ^ 2 ifc- 2 i'fc' ® ‘^ 2 iA:- 2 j'fc'’ (B.17), the range of indices spanned by 2^k — 2^'k' is 

Z gcd(2'^ , 2-^'). Thus, we would like to show that 

OO 

ll‘^ 2 gcd( 2 l, 2 J')ll < (B-21) 

Z= — OQ 

Since the matrix norm (see (2.2)) is sub-multiplicative, ll^ 2 gcd( 2 .:' 2 i') Hl = 

l|vec($^g^d(2^2i'))vec($^gcd(2i,2/))1bi ^ l|vec(^^g,d( 2 i, 2 i'))llp- Moreover, Proposition 3.2 

implies that the sequence { vec(<I>^g^jj. 2 j 2 i')) is summable (with > 2; see (2.6)); hence, 

(B.21) holds. The expression (3.11) is now a consequence of Lemma B.3. □ 

Proof of Theorem 3.1 The argument is reminiscent of those in Bardet (2000), pp.510-513, 
Bardet (2002), p.997, and Istas and Lang (1997), Lemma 2. For notational simplicity, we will 
restrict ourselves to the bivariate context (n = 2). The argument for general n can be worked out 
by a simple adaptation. 

The proof is by means of the Cramer-Wold device. Form the vector of wavelet coefficients 
F = (di (2^M), d2 (2^M),..., di (2^C K,,) , d2 (2^1, iL,-J;...; 

di (2-?''", 1), da(2^-, 1),..., di(2^™, ), d2(2'^"*, KjJ) G , 

where T(p) := 2 ^j- Notice that m, ji,... ,jm are fixed, but each Kj goes to infinity with 

p. Let 

= (B.22) 

where 

= («i,i,«i,i2,ai,2)* G j = 31, ■ ■ ■,3m- 

Now form the block-diagonal matrix 



We would like to establish the limiting distribution of the statistic 


= y ^vec5lT(2^') = Y*DY, 

^ V2J 
3=n 


31 



where it suffices to consider cx in (B.22) such that 


cx*T.{H,AA*)cx>{) (B.24) 

(see Brockwell and Davis (1991), pp. 211 and 214). The matrix 'S{H, AA*) in (B.24) is obtained 
from Proposition 3.3, and can be written in block form as 

= (B.25) 

corresponding to block entries of the vector a = (a^j,.. •, . Let P = Cov(y, Y) and consider 

the spectral decomposition P^/^DP^/^ = OKO*, where A is diagonal with real, and not necessarily 
positive, eigenvalues A and O is an orthogonal matrix. Now let Z ~ N{Q, Then, 

T{u) 

= Z*OAO*Z = Z*AZ =: ^ Xiii^)Zf. 

i=l 

Assume for the moment that 

max |Ai(z/)| = o(—* 7 ^). (B.26) 

By (B.24) and Proposition 3.3, 


Jm Jm 

i/Var(T.) = E E “l{ 


^ V ^Cov(vec5lP(2^),vec5 


Jm Jm 

1T(2^'))}«,, ^ E E > 0. 


Therefore, there exists a constant C > 0 such that, for large enough u, i^Ya,r{T^) > C > 0. In 
view of condition (B.26), 

maxi=i x(^) |Ai(z^)| ^ _ 1/2 

-^- < C u ' max |Aj(z/)| —)• 0, —)■ 00 . 

y'Var(Ti,) i=i,...,T{u) 

The claim (3.12) is now a consequence of Lemma B.4. 

So, we need to show (B.26). The first step is to establish the bound 


sup max —II^jIIp sup u^Tu. 




(B.27) 


Let u G ^ and let v = T^/^u. We can break up the vector v into two-dimensional subvectors 

u.,. to reexpress 

V = ..., ;...; ..., )*. 

Based on the block-diagonal structure of D expressed in (B.23), 


|^*pl/2^pl/2^| ^ |v*Dv| = lEEL 


" A'V2J 

1 = 11 1 = 1 J ' 


Jm . 

IK,I 

7 = 71 1 = 1 


j=jl,-,jm K 


where the constant C comes from a change of matrix norms (see (2.2)) and only depends on the 
fixed dimension n = 2. By taking sup^g_ 5 T(,/)-i on both sides of (B.28), we arrive at (B.27). 
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The second step towards showing (B.26) consists of analyzing the asymptotic behavior of the 
right-hand side of (B.27), as v ^ oo. Note that 


1 

max — 
Kj 


im-llh 


r\j 


C 


1 


V —)■ oo. 


(B.29) 


In view of (B.29), for establishing (B.26) it suffices to show that sup^g_ 5 T(i/)-i u^Tu is bounded. 
So, rewrite Y = (hi)j=i Since 


sup u*ru < 


T(^) 

■ V|Cov(yi,,yi2)| 

*1 = 1 ,... 


(B.30) 


(see Lemma 1 in Bardet (2000), p.509), we only need to show that the right-hand side of (B.30) 
is bounded. Fix 0 < d < 1/2 and without loss of generality assume that length(iF) = 1. We 
turn back to wavelet and dimensionality parameters (indices) to reexpress, and then bound, the 
right-hand side of (B.30) as 


im ^j' 2 


max max max 
k=l,...,Kj *=1,2 


EEE \Cov{di{2^ ,k'))\ 

j'=ji k'=l i'=l 




< 2m max max max 
k=l,...,Kj *,*'=1,2 


\CoY{di{2^,k),dA2^\k'))\ 

k'=l 


= 2m max max max 
k=l,...,Kj *,*'=1,2 




V1 f 

/ I max{21?,21? } 

k'=l \ |2ifc-2i'fc'| 




max{21? ,2J } 


<5 


|Cov(d,(2^fc),d,*(2^'',fc'))|. 


|2Jfe-2i fc'l “ J 


(B.31) 


One can interpret the bound in (B.31) as dividing up the covariance terms into those that, 
parameter-wise, are either far apart or close together (see (3.4)). 

We now develop a bound for the first summation term in (B.31). By the property {P2), 


\CoY{di{2^,k),di'{2^ ,y))\ < max max YYar di{2^,0)\/YaT dii{2^',0) < C, 

i.e., the covariance terms have a common bound. Moreover, in regard to the associated indicators 
in (B.31), when j' > j, #{/c' : \2^~^'k — k'\ < 5“^} < 2d~^ + 1. Alternatively, when j' < j, 
#{k' : \2^~^'k — k'\ < 2^~^'S~^} < 2 2^~^'+ 1. Therefore, the first summation term in (B.31) 
comprises finitely many terms and is bounded by a constant, irrespective of u. 

To bound the second summation term in (B.31), since \Coy{ di{2^,k),di'{2^ ,k'))\ < 

||Cov(Z)(2-^, A:), D(2'^/A:'))||p, the bound (3.6) implies that 


K,. 


El 


max{2-? ,2d 
k' = l 1 \23k-2i' k'\ ““ 


Lfc'l - J 


\CoY{di{2^ ,k),di,{2^' ,k'))\ 


<cj;i 

k'= 


log^ \2^k- 2^'k'\ 


/ f max{2i,2j'} <A \2j k - 2^ fe'I 

1 f \23k-2i'k'\ “ J 


2N^-2 - 




^^0 


(B.32) 
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where 2N^ — 2 > 1 by (2.6) and C does not depend on /c, k'. Consequently, (B.31) is bounded, 
and so is (B.30), as we wished to show. This establishes (B.26), and as a result, also (3.12). □ 

The next lemma, stated without proof, is used in the proof of Theorem 3.1. It is a simple 
consequence of the Lindeberg central limit theorem. 

Lemma B.4 Let {Wj^n}j=i,...,n, n G N, 6e an array ofi.i.d. random variables such that EWj^t = 0 
and EWjj^ < oo. Let {Xj,n}j=i,...,n, n G N, be an associated array of constants Xj^n £ 1^- Define 
the statistic Vn = ^j,nWj^n o,nd its variance = Var(14.)- |= o{an), 

then ^ A A^(0,1). 

C Estimation based on the discretized wavelet transform 

In this section, we address the issue raised in Remark 5.1. Instead of a continuous time OFBM 
path {Rff(t)}tgR, we assume that only a discrete OFBM sample 

{BH{k)]k&^ (C.l) 

is available. The main goal is to develop the asymptotic distribution of the Hurst matrix 
estimators starting from the so-called discretized wavelet coefficients (as defined in (C.5) below). 
The mathematical construction in this section is carried out along the same lines of that in Stoev 
et al. (2002) (notably Proposition 2.4 and Theorem 3.2). 

We suppose the wavelet approximation coefficients stem from Mallat’s pyramidal algorithm, 
under a multiresolution analysis of L^(M) (MRA; see Mallat (1999), chapter 7). Accordingly, we 
need to replace (1F2) with the following more restrictive condition. 

Assumption (1F2'); 

the functions (f (a bounded scaling function) and if correspond to a MRA of L^(M), 
and supp((/?) and supp(V’) are compact intervals. 

All through this section, we assume that (IFl), (1F2') and (1F3) hold. Given (C.l), we initialize 
the algorithm with the vector-valued sequence 

M” 9 ao,fc := a^Buik), /c G Z, := / ip{t)dt, 

Jr 

also called the approximation coefficients at scale 2^ = 1. At coarser scales Mallat’s algorithm 
is characterized by the iterative procedure 

^ ^ hy—2k^j,k' ^ d,jj^\ k ^ ^ 9k' — 2kOij^k' 1 j ^ ^ ^ 'L, 

k'ez fc'ez 

where the filter sequences {hk}k&z, {s'fcjfcez called low- and high-pass MRA filters, respec¬ 
tively. Due to (VF2'), only a finite number of filter terms is non-zero, which is convenient for 
computational purposes (see Daubechies (1992), chapter 6). 

For the subsequent developments, we will need to strengthen the condition (OFBMl). 

Assumption (OFBMl'): the eigenvalues hk of the matrix exponent H satisfy 
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G (0,l)\{l/2}, fc = l,...,n. (C.2) 

Let D be the shifted exponent matrix (1.3). Under (C.2), Theorem 3.2 in Didier and Pipiras 
(2011) yields the time domain representation 

= \ [ 9H{t,u)B{du)\ , (C.3) 

1 Jr J tm 

where B{du) is Brownian random measure such that EB{du)B{du)* = du and 

gH{t,u) = {{t - tt)+ - (-rt)+}M+ + {{t - u)^ - (-u)?}M_ G M(n,M) (C.4) 

for M+,M_ G M(n,M) (in particular, the theorem describes the relation between the matrix 
parameters M+, M_ and A in (1.2)). The next result draws upon (C.3) to provide an integral 
representation for the normalized discretized wavelet coefficients 

M*" 9 b{2^, k) := 2-B'^dj^k- (C.5) 

For notational simplicity, we define ipi{u) = ip{u — i), i G Z. 

Lemma C.l Fix j G N and k G Z. Under the assumptions (OFBMf) and (OFBM2), let 
D{2^,k) be the discretized wavelet coefficient (C.5). Then, 

p p OO . 

b{2\k)= / 2^^| / V’(i) X] -k'^dt^B{ds). (C.6) 

JM JM i=-oo 

Proof: The assumed conditions on ip and ip and the expression (C.4) imply that the integral 
on the right-hand side of (C.6) is well-defined in the mean-square sense. Therefore, (C.6) can be 
shown by a direct adaptation of the proof of Lemma 6.1 in Stoev et al. (2002), which also makes 
use of Lemma 6.2 in the latter reference. □ 

For the proofs of Proposition C.l and Theorem C.l below, recall the matrix norm notation 
(2.2). Then, for a random vector X, we can define the norm ||X||£,q(p) = (Fi||X||^,)^/'?, g G N. 

The next proposition establishes that the mean-square distance between the discretized and 
regular wavelet coefficients (expressions (C.5) and (1.5), respectively) is bounded by a constant 
that does not depend on the octave j or the shift k. 

Proposition C.l Under the assumptions (OFBMf) and (OFBM2), let D{2fk) he the approxi¬ 
mation coefficient (C.5). Then, there is a constant C{Hi, ip, (p) > 0, not depending on j ork, such 
that 

\\D{2^,k)-b{2^,k)\\L2^P)<C{H,iP,ip). 

Proof: Without loss of generality, we can assume that the interval [—K,K\, K > Q, contains 
the support of both functions ip and ^p. The relations (C.3), (2.6) and (A.l) yield the integral 
representation D{2fk) = ip{t)gH{t, ^ — k)df\B{ds). So, let 

„K °° 

/h(s) = 2^^ / ip{t) guipt,^ - k'^ - ^ a^ipi{2H)gH[^,-^ - l^ dtGM{n,R). 

i=—ao 
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Then, by expression (C.6), we can write 


= JtT{JjH{s)fH{syds] = J{JjfH{s)\\lds]. (C.7) 


However, for almost every s E 


pK ^ 

\\fHis)\\i 2 < 2^^ J - k^(l - a^(pi{2H)^dt 


rK oo 


+ 2^^ y{t) y] ^ -/c) ^ - A:))a<pV9i(2^t)dt . (C.8) 

i=-oo 12 

By Lemma 6.2 in Stoev et al. (2002), the first term on the right-hand side of (C.8) is zero. 
Therefore, by a change of variables and Fubini’s Theorem, (C.7) is bounded by 

C'|a^|2^/2( / 2^^ f yit) (^gH{t,s) - gH(^^,s'j'j(pi{2H)dt d.sj ^ . 


“ 2x1/2 


Y 


< E ||||^*-*(‘ 

L2(]R) i=-oo 


IfIIl2(]r) 


C|a^|2^'/2 


/ J\ ' 

^y{t){gH{t,s) - gH(^^,^^‘Pi{‘2.H)dt ds 


oo I.K 


<C\ay 2^^‘^ Y y - gH(^^,s^^ ^ds'j dt, (C.9) 

where the last inequality is a consequence of Lemma C.3 below. On the other hand, from (C.4), 
gH{t,s) — gui^ 1 s) = gH{t — ^,s — T). Thus, by making the change of variables u = s — if2y 


i i 


trf2^ gH[t- —,s- — ]gH[t- —,s- — 


21 21 




21 21 




= / tri2i gnit- —,u\gH(t- —,u 


^210*} 


du = tr< 2-^^ t-r S t — 


2^^*Y (C.IO) 


where the last step follows from the fact that EB}{{t)B}{it)* = 

|ty£'i?i^(sign(t))i?//(sign(t))*|ty* =: |tyS|ty*, t G M. By applying (C.IO) and account¬ 
ing for supp(V’) and supp((/?i(2'^t)), the expression (C.9) can be bounded from above by 


C|a^|2^'/2 


■^(22+1) 

E / 

i=-K{2l+l) h'i--K)/22 


tri 2 ^'^lt-yEt-y^^ 




= C|a^|2^'/2 


= C|a^|2^'/2 ii^ii^ ii^ii^ (2K(2^- + 1) + 1)2 f {tv{2^^2^^*}y/^dw 

Jo 

IWIloo llx'lloo l\tr{.«Ez«‘})'''‘dz < la^l ||^|U ll»»IL 


V\ IIV^Iloo ll'f'lloo 


C{H), 
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where the equalities result from the consecutive changes of variables w = t — ^, z = 2^ w. This 
proves the claim. □ 


The next definition describes the discretized counterparts of the sample wavelet variance (4.1) 
and the estimators (1.7) and (4.42). 

Definition C.l Let j, log 2 a(zz) G N, and let D{a{v)2^,k) be the discretized wavelet coefficient 
(C.5) at scale a{u)2K We define the associated sample wavelet variance as 


W{a{iy)2^) 


1 

— Y,Dia{,z)2\k)D{aiu)2^,kr. 
^“4 k=i 


Under (2.12), we also define the operator-renormalized discretized wavelet coefficients 

b^{2\k) := a{v)-^b{a{u)2^, k), D^{2\ k) := a{v)-^ D{a{v)2\k), (C.ll) 

and the random matrix 

Ka,j ^ 

B^{2b = P-^j^^D^{2^,k)D^{2^,ky{P*)-\ PeGL{2,R) (C.12) 

^“4 k=i 


(c.f. relation (4.10)). Let Ai(a(iz)2'^) < A2(a(jz)2-^) be, respectively, the smallest and largest 
eigenvalues of W(a(iy)2^). By analogy to Definitions 4.1 and 4.2, we define the estimators 


hi{a{u)2b 


logAi(a(;z)2-^') 
21oga(zz) ’ 


h2{a{v)2b 


log X 2 {a{iy) 2 b 
21oga(jz) ’ 


6{a{v)2b 

for hi, /i 2 and 9, respectively, where 


^i(«(^)2^) -QyaH 

bj,a{u) 


(C.13) 

(C.14) 


I yj,a{P 


:= W{a{v)2^) = Pdiag(a(zz)^^,a(i/)^^)i?iy( 2 -^)diag(a(zz)^^,a(jz)^ 2 ^p*_ (C.15) 


Theorem C.l below shows that, under mild conditions, estimating the Hurst eigenvalues based 
on either the discretized or the regular wavelet transform yields the same asymptotic behavior. 
Before stating and proving it, we will need the next lemma (c.f. Lemma 4.1). 


Lemma C.2 Under the assumptions (OFBMl'), (OFBM2) - (4) and (4.2), let By{2b, Ba{2^) 
be as in (C.12) and (4.10). Then, 


ybyyj{By2b - Ba{2b) 


LUP) 


—> 0 , 


oo. 


(C.16) 
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Proof: By Minkowsky’s inequality, the left-hand side of (C.16) is bounded by 


K„, 


C 


-j= V \\d,{ 2^, k)D,{2^, k)* - D,{2^, k)D,{2^, ky 

k=i 


LHP) 


(C.17) 


However, for k = 1,..., Ka,j, the deviation between the associated non-normalized wavelet vari¬ 
ance terms can be recast as 


D{2^a{u), k)D{2^a{u), kf - D{2^a{v), k)D{2^a{u), kf 

= [b{2^a{v),k) - D{2^a{v),k)][b{2^a{u),k) - D{2^a{v),k)]* 

+[b{2^a{v), k) - D{2^a{u), k)]D{2^a{u), k)* + D{2^a{u), k)[b{2^a{u), k) - D{2^a{u), k)]*. 

Therefore, by Proposition C.l, property (P3), and by using the fact that the 1 ^ matrix norm is 
sub-multiplicative, each term under the summation sign in (C.17) can be bounded by 

\\a{i')~^{[D{2^a{i'), k) — D{2^a{u), k)][D{2^a{v), k) — D{2^a{v), k)]* 

+ [b{2^a{u), k) - D{2^a{v), k)\D{2^a{v), k)* 

+D{2^a{v), k)[b{2^a{v), k) - D{2^a{u),k)]*}a{u)~^* WlHp) 

< II [b{2^aii^), k) - D{2^a{v), k)][b{2^a{v),k) - D{2^a{uy k)]* H^qp) 

+ ||a(z.)-^||p||5(2Mz.),fe)-Zl(2%(^.),fe)||pi(p)||D(2^fe)||pi(p) 

+ ||75(2^fc)||ii(P)||5(2Mi^),fc)-Zl(2M^),A:)||ii(p)||a(z.)-^|bi 

< l|a(i^)“^||p||7)(2%(i/),/c) -Il(2^a(z/),/c)|||2(p) 

+ ||a(z/)-^||p||5(2M^^),A:)-Zl(2M^^),fc)||p2(p)||Zl(2^0)||ii(p) 

+ ||Zl(2^0)||ii(p)||5(2Mi^),A:)-71(2^a(z/),A:)||p2(p)||a(i/)-^*||p 

< ||a(i/)-^||2 c2(77, y, + ||a(^.)-^||p ||Z1(2^ 0)||ii(p) C{H, y, 

+ ||a(z.)-^lp||D(2^0)||ii(P)|| C{H,y,ip) < C||a(z.)-^||p. 

Consequently, (C.17) is bounded by C'Y^ifa,j||a(z^)~'^||p —)• 0, i/ —)• 0, where the null limit results 
from (4.2). This establishes (C.16). □ 


Theorem C.l For m ^ N, let ji < ... < jm be a set of fixed octaves j. Let hi, h 2 , hf, be as 
in (C.13) and (4.9). Under the assumptions (OFBMf), (OFBM2) - (4) and (4.2), asv^oo, 


(^) 


2\og{a{v)2L)^Ka,j[hi{a{v)2b - hf{a{v)2b\ 
2\og{a{v)2L)^/K~[h2{a{v)2b - h^{a{v)2b] 


—)■ N (0, ^hi,h2 (il) • • • ) jm)), 


(C.18) 


where the matrix T,h^^h2{ji^ ■ ■ ■ , jm) is given in Theorem 4-1; 


(ii) if, in addition, P 22 0 and the condition (4.45) holds, then 


a{uf^-^^ybybj{0{a{n)2L) - 6{a{u)2b} 4 N{0,al), n ^ 00 . 
where the variance is given in Theorem 4-2. 


(C.19) 
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Proof: We begin by showing (i). Rewrite the left-hand side of (C.18) as 


f \/^<^(logAi(a(jv) 2 ^) - logAi(a(jv) 2 ^))-h y^^(logAi(a(i/) 2 ^) - logAf(a(z^) 2 ^)) \ /q 20 i 
V \/^^(^osM(a(iy)2^) -logA 2 (a(iy) 2 ^)) + ^/7^(log^2(a(ly)2J) -logAf(a(iy)2^)) j . 

where j = ji, ■ ■ ■, jm- As in the proof of Theorem 4.1, without loss of generality we can fix j. 
Theorem 4.1 implies the asymptotic normality of the second term in the sum (C.20). We will 
show that the first term in the sum (C.20) goes to zero in L^{P). Note that 

2log{a{v)2^)\hi{a{v)2^) -hi{a{iy)2^)] = logXi{a{u)2^) - log \i{a{iy)2^), i = 1,2. (C.21) 


So, we can repeat the proof of Theorem 4.1 with W{a{u)2^) in place of EW{a{v)2^) and W{a{v)2^) 
(see (C.15)) in place of W{a{v)2X). In regard to i = 1 in (C.21), we arrive at a relation analogous 
to (4.24), where 


_ ac — b'^p 
r := 4--^ ~ 


4(det(P))2det By{2^) 


(a+ 2)2 ((^ 2 ^ 

The analogous expression to (4.29) is 


V —)■ oo. 


- r (ac — 6^) 

ni 

1 

to 

(det(P))2 f 


(a + c)2 

(a + c)2 J 

a(p)2(^2-fti) 1 


r 


Op(l) 


Op(l) 


0 . 


\Op{a{uY^ ^^) + {pl2+pl2)b22f[Op{a{u)^^ '*") + (p?2 + ^ 22 )^ 22 ]^^ ^ 

(C.22) 

The numerators of the first and second terms appearing in the sum in (C.22) are bounded in 
probability as a consequence of Lemma C.2 and of applying the mean value theorem to the 
function fz{x) = x^. The relation analogous to (4.31) is 


detR^(2^) - det.Ba(2^) 622 - 622 + ^ 

det.Ba( 2 -?') 622 

where the limit follows from Lemma C.2. In regard to i = 2 in (C.21), it suffices to repeat the 
argument in Theorem 4.1 to arrive at the relation 

V^{log(A 2 (a(p) 2 ^)) - log(A 2 (a(zz) 2 ^))} ^ ^ 22 ( 2 ^) - 622 ( 2 ^) \ P 

^ 022 ( 2 -^) 2 

where the limit is again a consequence of Lemma C.2. Therefore, (C.18) holds. 



We now turn to {ii). Recast the left-hand side of (C.19) as 

- 0{a{yW)}. (C.23) 

The weak limit of the second term in the sum (C.23) is given by Theorem 4.2. By analogy with 
(i), we will show that the first term in (C.23) converges to zero in probability. Rewrite this term 
as 




h 


^\+a{up-^^ 




(C.24) 
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(C.25) 


~ • ■ P 

As a consequence of the proof of (C.18), log(Ai(a(z/)2'^)/Ai(a(i/)2'^)) —)■ 0, namely, 

Ai(a(i/)2'^)/Ai(a(i/)2'^) 1. 

Based on (C.25), Lemma C.2 and the expression (C.15), we obtain the relations 

Ai p 1 a{vY^^2{det{P)Y det By{2^) 

b [Op(a(p)^i-'^2) + pi2P22b22\a{i^y^^ [Op(a(p)'^i-'^2) + (^ 12 +^ 22 )^ 22 ]’ 

| = Op(l), 7^1^} ^0. (C.26) 

Therefore, 

a(O'‘»-'«(-l))^N/5f^{^}^0. (C.27) 

0 ^ b 

Moreover, for y/Kaj{Xi — Ai}, we can obtain relations analogous to (4.51), (4.52), (4.53) and 
(4.54) by making use of Lemma C.2. We conclude that 

To show that _ 

a(O^-+yiL-{|- 2 }^ 0 , 

'' b b ■’ 

we can follow the steps of the proof of the analogous statement in Theorem 4.2 to arrive at 
expressions which, in turn, are analogous to (4.59) and (4.60). Then, we can make use of Lemma 
C.2 to find a zero limit in probability, instead of a distribution. Thus, (C.19) holds. □ 


The next lemma, a multivariate generalized Cauchy-Schwarz inequality, is used in the proof of 
Proposition C.l. It can be shown by an adaptation of the univariate argument, which we provide 
for the reader’s convenience. 


Lemma C.3 Let g : 


Then, 


M(n,M), {t,s) I—)■ g{t,s), be a function such that 




g{t, ■)dt 






g{t,s)dt 


ds 


1/2 


< 


|| 5 ((t,s)||p ds) ^ dt. 


(C.28) 

(C.29) 


Proof: Starting from the square of the left-hand side of (C.29), 


g{ti, s)dti 


g{t2,s)dt2 


ds < 




g{ti,s)dti 


P 


)(^ 115(^2, s)||;2dt2)ds 


< 




\\git 2 ,s)\\p ds]dt 2 


dsi) ^ J 115(^2, S 2 )||p dS 2 ) ^ dt 2 . 


The equality is a consequence of Fubini when applied to the integrators dt 2 ds, whereas the second 
inequality is a consequence of (C.28) and the Cauchy-Schwarz inequality. This establishes the 
claim (C.29). □ 
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